{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "# 3.3 Discontinuous Galerkin Discretizations\n",
    "We are solving the scalar linear transport problem\n",
    "\n",
    "Find $u: [0,T] \\to V_D := \\{ u \\in L^2(\\Omega), b \\cdot \\nabla u \\in L^2(\\Omega), u|_{\\Gamma_{in}} = u_D\\}$, s.t.\n",
    "\\begin{equation}\n",
    "\\int_{\\Omega} \\partial_t u v +  b \\cdot \\nabla u v = \\int_{\\Omega} f v \\qquad \\forall v \\in V_0 = \\{ u \\in L^2(\\Omega), b \\cdot \\nabla u \\in L^2(\\Omega), u|_{\\Gamma_{in}} = 0\\}\n",
    "\\end{equation}"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "outputs": [],
   "source": [
    "import netgen.gui\n",
    "%gui tk\n",
    "from math import pi\n",
    "from ngsolve import *"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "As a first example, we consider the unit square $(0,1)^2$ and the advection velocity $b = (1,2)$. Accordingly the inflow boundary is $\\Gamma_{in} = \\{ x \\cdot y = 0\\}$. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {
    "slideshow": {
     "slide_type": "-"
    }
   },
   "outputs": [],
   "source": [
    "from netgen.geom2d import SplineGeometry\n",
    "geo = SplineGeometry()\n",
    "geo.AddRectangle( (0, 0), (1, 1), \n",
    "                 bcs = (\"bottom\", \"right\", \"top\", \"left\"))\n",
    "mesh = Mesh( geo.GenerateMesh(maxh=0.2))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "We consider an Upwind DG discretization (in space):\n",
    "\n",
    "Find $u: [0,T] \\to V_h := \\bigoplus_{T\\in\\mathcal{T}_h} \\mathcal{P}^k(T)$ so that\n",
    "\n",
    "$$\n",
    "  \\sum_{T} \\int_T \\partial_t u v + b \\cdot \\nabla u v + \\int_{\\partial T} b_n (\\hat{u}-u) v = \\int_{\\Omega} f v, \\quad \\forall v \\in V_h.\n",
    "$$\n",
    "\n",
    "Here $\\hat{u}$ is the Upwind flux, i.e. $\\hat{u} = u$ on the outflow boundary part $\\partial T_{out} = \\{ x\\in \\partial T \\mid b(x) \\cdot n_T(x) \\geq 0 \\}$ of $T$ and $\\hat{u} = u^{other}$ else, with $u^{other}$ the value from the neighboring element."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "-"
    }
   },
   "source": [
    "There is quite a difference in the computational costs (compared to a standard DG formulation) depending on the question if the solution of linear systems is involved or only operator evaluations (explicit method). We treat both cases separately:\n",
    "\n",
    "* Case 1: Solution of the time-dependent problem by explicit time stepping\n",
    "* Case 2: Solution of the stationary problem "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## Case 1: Explicit time stepping with a DG formulation\n",
    "\n",
    "Explicit Euler:\n",
    "\\begin{align}\n",
    "  \\sum_{T} \\int_T u^{n+1} v & = \\sum_{T} \\int_T u^{n} v - \\Delta t \\sum_{T} \\left\\{ \\int_T  b \\cdot \\nabla u v + \\int_{\\partial T} b_n (\\hat{u}-u) v \\right\\} - \\Delta t \\int_{\\Omega} f v, \\quad \\forall v \\in V_h, \\\\\n",
    "  M u^{n+1} &= M u^{n} - \\Delta t C u^n + \\Delta t \\ f \n",
    "\\end{align}\n",
    "\n",
    "In our first example we set $u_0 = f = 0$. "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "### Computing convection applications $C u^n$\n",
    "* We can define the bilinear form **without setting up a matrix** with storage.\n",
    "* A `BilinearForm` is allowed to be **nonlinear in the 1st argument** (for operator applications)\n",
    "* For the DG formulation we require integrals on element boundaries, keyword: **`element_boundary`**\n",
    "* To distinguish inflow from outflow we evaluate $b \\cdot \\mathbf{n}$. Here the normal $\\mathbf{n}$ is available as a `specialcf`\n",
    "* To make cases with `CoefficientFunctions` we use the `IfPos`-`CoefficientFunction`. **`IfPos(a,b,c)`**. `a` decides on the evaluation. If `a` is positive `b` is evaluated, otherwise `c`.\n",
    "* To access the **neighbor** function we can use `u.Other()`\n",
    "* To incorporate **boundary conditions** ($\\hat{u}$ on inflow boundaries), we can use the argument `bnd` of `.Other(bnd)`. If there is no neighbor element (boundary!) the `CoefficientFunction` `bnd` is evaluated. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "outputs": [],
   "source": [
    "b = CoefficientFunction((1,2))\n",
    "n = specialcf.normal(mesh.dim)\n",
    "ubnd = IfPos(x,1,0)\n",
    "\n",
    "V = L2(mesh,order=2)\n",
    "u,v = V.TrialFunction(), V.TestFunction()\n",
    "\n",
    "c = BilinearForm(V)\n",
    "c += SymbolicBFI( b * grad(u) * v)\n",
    "c += SymbolicBFI( IfPos( (b*n), 0, (b*n) * (u.Other(ubnd)-u)) * v, element_boundary=True)\n",
    "\n",
    "gfu_expl = GridFunction(V)\n",
    "Draw(gfu_expl,mesh,\"u_explicit\")\n",
    "\n",
    "res = gfu_expl.vec.CreateVector()\n",
    "c.Apply(gfu_expl.vec,res)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "### Solving mass matrix problems\n",
    "\n",
    "* Need to invert the mass matrix\n",
    "* For DG methods the mass matrix is block diagonal, often even diagonal.\n",
    "* In `NGSolve`, `FESpace`s can provide a routine `SolveM` which realizes the operation $M^{-1} u$. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {
    "slideshow": {
     "slide_type": "-"
    }
   },
   "outputs": [],
   "source": [
    "t = 0\n",
    "dt = 0.001\n",
    "tend = 1\n",
    "\n",
    "while t < tend-0.5*dt:\n",
    "    c.Apply(gfu_expl.vec,res)\n",
    "    V.SolveM(res,CoefficientFunction(1.0))\n",
    "    gfu_expl.vec.data -= dt * res\n",
    "    t += dt\n",
    "    Redraw(blocking=True)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## Case 2: Solving linear systems with a DG formulation\n",
    "\n",
    "When it comes to solving linear systems with a DG formulation we have to change the way the sparsity pattern is typically constructed."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "### Sparsity patterns in NGSolve\n",
    "* sparsity pattern of a standard FEM and DG matrix is different\n",
    "* In `NGSolve` the sparsity pattern is set up whenever a `BilinearForm` is assembled\n",
    "* the sparsity pattern is based on the finite element space only (and not the integrators!).\n",
    "\n",
    "Two cases:\n",
    " * standard FEM formulation (element-based couplings only)\n",
    " * DG formulations (element- and facet-based couplings) (`dgjumps=True`)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "#### 1. Standard formulation\n",
    "* In a standard formulation unknowns only have a nonzero entry if the support of corresponding basis functions overlap\n",
    "* In `NGSolve` this sparsity pattern is obtained by allocating nonzero entries whenever two unknowns have an association with the same element.\n",
    "\n",
    "The idea is sketched here for one dof:\n",
    "\n",
    "dof $\\to$ elements of dof = [el$_1,\\dots,$el$_n$] $\\to$ dofs of (el$_1$) $\\cup \\dots \\cup$ dofs of (el$_n$)."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "#### 2. DG formulation ('dgjumps')\n",
    "* In a DG formulation couplings (nonzero entries) are introduced also for basis functions which do not have overlapping support.\n",
    "* An additional mechanism (keyword: `dgjumps`) has been introduced to determine the sparsity pattern. \n",
    "* This can be activated by adding the `dgjumps` flags to the `FESpace`. \n",
    "* We do not need the flag if we do not set up matrices (explicit time ste\n",
    "pping).\n",
    "* A `dgjumps`-formulation introduced additional couplings through facets, i.e. for all dofs that associate to the same facet a nonzero entry is reserved.\n",
    "\n",
    "The idea is sketched here for one dof:\n",
    "\n",
    "dof $\\to$ facets of dof = [fac$_1,\\dots,$fac$_n$] $\\to$ dofs of (fac$_1$) $\\cup \\dots \\cup$ dofs of (fac$_n$).\n",
    "\n",
    "Note that 'dofs of facet $F$' is always larger than 'dofs of element $T$' if $F \\subset \\partial T$."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "We want to demonstrate the difference in the sparsity pattern with a simple comparison in the next block"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {
    "slideshow": {
     "slide_type": "-"
    }
   },
   "outputs": [],
   "source": [
    "try:\n",
    "    V1 = L2(mesh,order=0)\n",
    "    a1 = BilinearForm(V1)\n",
    "    a1.Assemble()\n",
    "\n",
    "    V2 = L2(mesh,order=0, dgjumps=True)\n",
    "    a2 = BilinearForm(V2)\n",
    "    a2.Assemble()\n",
    "\n",
    "    import scipy.sparse as sp\n",
    "    import matplotlib.pylab as plt\n",
    "    plt.rcParams['figure.figsize'] = (12, 12)\n",
    "\n",
    "    a1.mat.AsVector()[:] = 1 # set every entry to 1\n",
    "    rows,cols,vals = a1.mat.COO()\n",
    "    A1 = sp.csr_matrix((vals,(rows,cols)))\n",
    "    a2.mat.AsVector()[:] = 1 # set every entry to 1\n",
    "    rows,cols,vals = a2.mat.COO()\n",
    "    A2 = sp.csr_matrix((vals,(rows,cols)))\n",
    "    fig = plt.figure()\n",
    "    ax1 = fig.add_subplot(121)\n",
    "    ax2 = fig.add_subplot(122)\n",
    "    ax1.spy(A1)\n",
    "    ax2.spy(A2)\n",
    "except ImportError:\n",
    "    pass"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "#### non-`dgjumps` vs `dgjumps`"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {
    "slideshow": {
     "slide_type": "-"
    }
   },
   "outputs": [],
   "source": [
    "try:\n",
    "    plt.show()\n",
    "except ImportError:\n",
    "    pass    "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "* when assembling a `BilinearForm` the bilinear form has to be **linear** in both arguments again (Linearizations will be addressed later)\n",
    "* `...Other(bnd=..)` does not make sense any more\n",
    "* when setting up the system we require `LinearForm` and `BilinearForm`s separately to implement boundary conditions."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {
    "slideshow": {
     "slide_type": "-"
    }
   },
   "outputs": [],
   "source": [
    "b = CoefficientFunction((1,2))\n",
    "n = specialcf.normal(mesh.dim)\n",
    "ubnd = IfPos(x,1,0)\n",
    "\n",
    "V = L2(mesh,order=2, dgjumps=True)\n",
    "a = BilinearForm(V)\n",
    "u,v = V.TrialFunction(), V.TestFunction()\n",
    "a += SymbolicBFI( b * grad(u) * v)\n",
    "a += SymbolicBFI( IfPos( (b*n), 0, (b*n) * (u.Other()-u)) * v, element_boundary=True)\n",
    "a.Assemble()\n",
    "\n",
    "f = LinearForm(V)\n",
    "f += SymbolicLFI( (b*n) * IfPos( (b*n), 0, -ubnd) * v, BND, skeleton=True)\n",
    "f.Assemble()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "outputs": [],
   "source": [
    "gfu_impl = GridFunction(V)\n",
    "gfu_impl.vec.data = a.mat.Inverse() * f.vec\n",
    "Draw(gfu_impl,mesh,\"u_implicit\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "### Remarks on sparsity pattern in NGSolve"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "source": [
    "#### Remark 1: The sparsity pattern is set up a-priorily\n",
    "* The sparsity pattern of a sparse matrix in NGSolve is independent of its entries (it's set up a-priorily). \n",
    "* We can have \"nonzero\" entries that have the value\n",
    "\n",
    "Below we show the reserved memory for the sparse matrix and the (numerically) non-zero entries in this sparse matrix. \n",
    "\n",
    "(you need matplotlib/scipy for this block - but can ignore it otherwise)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAskAAAFdCAYAAADxKvAlAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMi4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvhp/UCwAAIABJREFUeJzt3Xv4XVV95/HPl5DESzSIAUwDTGKJY/GG8nusVu1QLVbBaXRqAdsqOoyxLdZa6zONlhmpjTOh1draCzYWClgr4p0SKkbEWq2oiXJP0dTEB9IIxEsEHU0D3/ljr0NOds75nX322Ze19n6/nidPztnn8ltn77W/53v2/q61zd0FAAAA4IDD2m4AAAAAEBuSZAAAACCHJBkAAADIIUkGAAAAckiSAQAAgBySZAAAACCn9STZzF5gZreb2XYzW9d2e8Yxs51mdrOZ3WBmW8KyI81ss5l9Pfz/qAjaebGZ3W1mtwwtG9lOy7wrrPubzOxpEbX5fDPbFdb3DWZ22tBjbwptvt3MfqGlNh9nZteZ2W1mdquZ/XZYHu26nqfNsa/rh5jZl8zsxtDuPwjLV5nZF0P7PmBmi8LyxeH+9vD4yjba3VWpxGwpjbhNzG6szcnF7AntjnZ9dypmu3tr/yQtkPRvkh4raZGkGyWd2Gab5mnrTknLcsv+SNK6cHudpAsiaOfPSnqapFsmtVPSaZL+UZJJeoakL0bU5vMlvXHEc08M/WSxpFWh/yxooc3LJT0t3H6EpK+FtkW7rudpc+zr2iQtCbcXSvpiWIdXSDorLH+3pN8It39T0rvD7bMkfaDpNnf1X0oxO7Q3+rhNzG6szcnF7AntjnZ9dylmt30k+emStrv7N9x9n6TLJa1puU3TWCPp0nD7UkkvbrEtkiR3/6yk7+QWj2vnGkmXeeZ6SUeY2fJmWnrAmDaPs0bS5e7+Y3ffIWm7sn7UKHff7e5fCbfvlbRN0gpFvK7nafM4saxrd/f7wt2F4Z9Leq6kD4Xl+XU92AYfkvQ8M7OGmtt1qcdsKbK4TcxuRooxO7Q1ubjdpZjddpK8QtIdQ/fv1Pwbv00u6ZNmttXM1oZlx7j77nD7W5KOaadpE41rZ+zr/7XhNNfFQ6dEo2tzODX0VGW/lpNY17k2S5GvazNbYGY3SLpb0mZlR0e+5+77R7TtwXaHx/dKenSzLe6saPpEQanG7STiyAhRx5GBFGO2lFbc7krMbjtJTsmz3f1pkl4o6Vwz+9nhBz07TxD9Nb5TaaekCyX9pKSTJO2W9I52mzOamS2R9GFJr3f37w8/Fuu6HtHm6Ne1u9/v7idJOlbZUZHHt9wkpCH5uJ1CG4Po44iUZsyW0ovbXYnZbSfJuyQdN3T/2LAsOu6+K/x/t6SPKtvodw1Ov4T/726vhfMa185o17+73xV2sgckvUcHThdF02YzW6gsaL3P3T8SFke9rke1OYV1PeDu35N0naRnKjv9eXh4aLhtD7Y7PL5U0rcbbmpXRdcn5pNw3I46joySQhxJMWZLacft1GN220nylyWtDiMeFykr2L6y5TYdwswebmaPGNyW9HxJtyhr69nhaWdL+ng7LZxoXDuvlPSKMIr3GZL2Dp12alWu9uslyta3lLX5rDAadpWk1ZK+1EL7TNJFkra5+58MPRTtuh7X5gTW9VFmdkS4/VBJpyqry7tO0kvD0/LrerANXirp0+EIEWaXRMyWko/b0caRcRKII8nFbCnNuN2pmJ0fydf0P2UjSL+mrF7l99tuz5g2PlbZaNEbJd06aKeymplrJX1d0qckHRlBW9+v7NTLfyir+TlnXDuVjUD9y7Dub5Y0F1Gb3xvadJOyHWj50PN/P7T5dkkvbKnNz1Z2Wu4mSTeEf6fFvK7naXPs6/rJkr4a2neLpP8dlj9WWfDfLumDkhaH5Q8J97eHxx/bRru7+i+FmD3UP6KP28TsxtqcXMye0O5o13eXYraFBgIAAAAI2i63AAAAAKJDkgwAAADkkCQDAAAAOSTJAAAAQE40SfLQ1ZCSQZubk2K7aXMzUmxzF6S43mlzc1JsN21uTirtri1JNrMXmNntZrbdzNYVeEkSKyyHNjcnxXbT5mak2OboELOjlWKbpTTbTZubk0S7a0mSzWyBsvkFXyjpREkvM7MT6/hbAIDZELMB4FC1zJNsZs+UdL67/0K4/yZJcvf/O+r5Cx621A9fenSh937SiqWF23Hzrr2FnzvN+5Z5/ybN91nuueceHXXUUQ22phoptLtsfyjT9+qSwnrOi6HNW7du3ePuaa24IXXGbEk6/DDTTy1/5MTnzRpTJ+1L99xzj761b9FMf6MNRddfTGLYLyeZpb8N97Vtu7+v/Q9Uk0tN+32Qwnoepe12F43Zh096QkkrJN0xdP9OST89/IRQj7JWkhY95gQtP/tPC73xHknLlizSlvNOnfjcles2FWyutGXD6YWfW+b9m7RnnsdMkgquP0ynbH8o0/cQFzP7ZtttmFFtMftBBeLOrDG1yL4Ua9yeZDium6QdxI2ZzdIXhvtalX2K74NmFI3ZdSXJE7n7RkkbJWnx8tVT/QTbc9++WtrUF6w/ANOaJWZLzcSdQbJS9EBKqrhObn/Mrd889b7T9f7fpLqS5F2Sjhu6f2xYVpm59ZvpBDMY/uXLDgX0Xu0xW8riThPxpg8HAoocvSS2p69MX+5D/29KXUnylyWtNrNVygLtWZJ+pco/sOe+fY0F3K4brMv5sJ6BTqs9Zg8Qu5tDbAdmU0uS7O77zey1kq6RtEDSxe5+ax1/q6pfTPMFEoIIv0y7oMxpO4n+3wdNxuwBYkoc2A7AeLXVJLv71ZKuruv9h9V9VGJcEFm2ZFGvAgxHJNJWtq/2qY/3WZMxe6CNI8p9i9tFENuB0VobuFe1UUGvymBY9EhzqiOnq9D3L54y/W3ZkvSmo6oTR7v7Z3h7E7Pj1NXYXra/EbcP6HrM7kySLB06mK/u6YYGuhpAypi0Trs8dVEKO3zsONrdb1vOO7X0l+406C+QDo3ZRftekVrvvuh6zO5Ukrznvn2HJMpNBNxhnMqbn+vQRDqVX5TTavIXdtP9HKhL0/2YmF1MH2ZEiqEfkHzHpVNJsnTwL7w2gl/ZX6Z9VvZXeeyBepZf2MPro8jRd/oYUjYcs5s2KoYQt+fX1ZgN5HUuSR4WQ5Cr+sqAOKAv0xtx4QD0RQwxWyJu16UvMXtg54SDG/Sh+HU6SU4Fp/vqw3ptRtcHb+Bghx9mbTehdcTterBOm0HMLiaKJPlJK5ZKPQ44RTvcqnWbOKpYQp+OXLSl64M3cLCfWv7IXsdsqVjcpmyjHGJ2/YjZxUSRJEtKfjqece2ucmcuMitEquuvTX3b6TFe2aRm4dGrnlJDc6JGzJ6Mso16ELMxUHfMjiZJHujaKSxGaqehDyO3y+rTnKBl9x07bEF0sbQpXYs5xOw0TDu4uU+I2ZMVjdnRBfYt553KL+sZ5JM71uX0+vaFlR9cki/rYU5QzIeYPRti9uz6VoaYj9n5o6nE7OpElyRL/LKuEuuynD7XxLX5hcMVsNJEnKkO67KcPl/Iqs3+0vWYHWWSnHqtW0yoiasHv9Tr0dUfHl1HzK4OMbseoy5khdl1PWZHmSQP8Iu6GWXXs6l/p7naMukqhV3/NY80ELObwXpOw6SxLmW2IzG7WVEnyVy9rhl1/hLs+zabVDtWVv49uv5rHmmgvrYZde7vbLOD43ZdMVtKN25POmjTJVEnyXkMEElP3780+/Z5ASB1xO3pdPlAWBRJ8s279o7slFWdnmhbn351TZLi9kN9OHWYpnExWzo4vqW6vxdJkroax1PdZkAdokiSx9lz377oTtcPTsPM8kszps9TVFWXsGRQSnflS0vQTzHG7TrE/vn6sA2AukWdJEvxBqK+/dpu8hKWfVu3deFLcnpl+54/cP/+GpqTLPpd+9gG6SFmT6/umB19khwr6qPrw9HmahBsp1f29Lld8KIbK24KgJ5JOWa3ddXaumM2SfIUmqyfrKq8oas42hzXZVmpLUbfEbMxSV8O7nTpu5kkuYA6ai3r2lm61DnnU/RLJX+J5a6qeqJ86ouByarY57oUsyfFDcoJkBqS5ALYsdOVP7ral1/yXcYROyBNlNL1V1UHrJqO4yTJBZAgdwdlGvVoMnFtchApgGYRo7upqjO6TfeNTiXJOzeczlFfTf4l3ucjakU+N31o+kEYJK5AeZNitqm6JCN2XGm3PsPlMByxL6ZTSbJ0YAejA4xHwJkfpwQPtue+fVwQB2hRXxLkUYjHaerKNulckgygfqMSZwCYBuMLmkXMnl5nkuRlSxZFcVqmyKwAKXbUsnViKUwNFkO/6bO25tcEBibF7RRjdgpmKdNim6AJySbJo4JaDDvNynWbOvlF37XPM4wEOR51bwvq9TFKl3+oNVWHmsIBEXRD/hoBZcqRFh696ilFnpdskhyzwa/crgVbzC7fJziKHZciR6jYr7uN/XG04WR7VNxi30Ebytbr22ELCuW/JMk1Itj2W5EvhFGPx3BGpE2xf372a/Rdl2ezycft2OMR6kWS3IJZ6nvZedNRZBtzJBmIHzG7P4jHGBZFkvykFUu1ZcPpvQkgnG7CAAEZKSJmp6OKbUS9MfoqiiQZAAC0o8isTH3FFQD7LaokuWhnTOlX7dz6zbUdhejytGyoBldYqkd+XS56zAknt9SUVnU1gagrbpddX3264l5sGDfSb1ElySmf0hqnzi+QLq4vVItgjjp1NYGoK27XGbO7sN6B2ESVJAOoH0eXAcSgq2ci0B0zJclmtlPSvZLul7Tf3efM7EhJH5C0UtJOSWe4+3dna2Yx7HAog36DPoktbqO/yh5ZZ1YgNKWKI8k/5+57hu6vk3Stu28ws3Xh/u9V8HcmGuxwHB3DNIoEavoUOiaauA1Mi5iNptRRbrFG0inh9qWSPiOCLQDEjLiNwmY5+8ZV+ZCSWZNkl/RJM3NJf+3uGyUd4+67w+PfknTMqBea2VpJayXp+OOPn7EZB+P0OapWpk8xiwgiVSpu1xmzkZY6j+RW9d3N7E+owqxJ8rPdfZeZHS1ps5n96/CD7u4hEB8iBOaNkrR4+WoftUOV/UU56jVt1jDlPxu/lNMT2zR+wAxKxe2iMbsrhj8fMTs9RbbXqLxgz337KNXAg2ZKkt19V/j/bjP7qKSnS7rLzJa7+24zWy7p7rLvX2XyEFO9MkkRBsp+8cbQj5GmOuP2nvv2dfKHXwqfhyOn05tlu057ARZidppKJ8lm9nBJh7n7veH28yW9VdKVks6WtCH8//EqGgqgek0nNPkvFr44mtVE3J7vhx/buz4c6QaqN8uR5GMkfdTMBu/z9+7+CTP7sqQrzOwcSd+UdMbszQRQB45k9w5xG0hY2dIfYnY5pZNkd/+GpKeMWP5tSc+bpVHDZrk8aMxzKVKn3D1l+lvZ7c6pVZTRRNwexLauxbQySUbX1gHKmyVm04faE/0V92ZJcqd57c4Np7f6SyvWZB7FldmGZbd7V4JmF+tXkal7u7Yds4ugb2OgKzG7b6JPkgF0V9kR6ACAdkwatNilmE2SDCBqVR2Bif2oIwDErkjJ3qiYnWr8JUkewqlfoBguroIyqo6xxOx4MW4iLmW2h0naMeVUd12TRJI8y+C9abR96pfLdSIVKfZDEqrm5WP3pH4z7dGmce8X0+nevsb1Ln6mlLW9PUzZpT7rfk3VkkiSYwl20sEdrenTBzGtByA1dX5J2AUv2lrbmyesrZiV0lzNxHX0QWxHpIvG7MPqbggAAACQGpJkAAAAICeJcotYUeMIAO0oe0Em4jaa0OTFpVAfkuQZdGmaEwBIWdGEJB+3idmoQ5MXl0J9KLcAAAAAckiSAQAAgBySZAAAACAnmZrkles2TV3UXnSARkpX+OnrxPQp4Cp0wKGmid19HVRHXG8HVwXEJMkkydL0Re19DCp9/IKJRR/7G1BE2UF1EgPrJOJ6XYjZmIRyCwAAACCHJBkAAADISarcAgDQH3PrN7fdBGCiMhcOQRo4kgwAiBKJB1JAP+0ukmQAAAAghyQZAAAAyCFJBgAAAHKSG7hX5qIiTYphMnwmpgcQm7n1m2uNO2Uv8BBDzC5iOK6bpB0bTm+vMagFFymJT3JJshR3kXz+SyDGifBjXn8AuqnuuFM2AU/xAibedgNQmZ382Ika5RYAAABADkkyAAAAkJNkuQUAAMOGSyQYd4GuKHOhEvp/dTiSDADoFMZdoCvK9GX6f3VIkgEAAIAckmQAAAAghyQZAAAAyEl24F7dE9NXJdaJ6hnkAqBp014Mqq34GWvcHlZkLmdiezPK9hcuHhK/ZJPk2APYQJEA1fbE9amsSwDpmybetHVxphTidhHE9mbwQ6S7KLcAAAAAckiSAQAAgJxkyy2k6evbMN6kU4esZwBVIXY3h9gOlDfxSLKZXWxmd5vZLUPLjjSzzWb29fD/o8JyM7N3mdl2M7vJzJ5WZ+Mlaq6awnoG0hF73JaIKbFgOwDjFSm3uETSC3LL1km61t1XS7o23JekF0paHf6tlXRhNc0EAEzhEhG3AWAmE5Nkd/+spO/kFq+RdGm4famkFw8tv8wz10s6wsyWV9VYAMBkxG0AmF3ZgXvHuPvucPtbko4Jt1dIumPoeXeGZYcws7VmtsXMtpRsAwCguJniNjEbQN/MPHDP3d3MvMTrNkraKEmLl6+e+vVVmFu/uVA9Vt0DG1KYuB7pK9rf8xjY0z1l4vZwzH74ise1ErNNUpE/bHU3RMRtNKNon8+/BtUomyTfZWbL3X13OC13d1i+S9JxQ887Niyb15NWLNWekg2ZRdEAV3cgLJqApDB5PeJVth+TCHRGZXH7p5Y/spWYXTRZaCKD78oFRxC3Mn25lV+wHVW23OJKSWeH22dL+vjQ8leE0dLPkLR36PQeAKA9xG0AmMLEI8lm9n5Jp0haZmZ3SnqLpA2SrjCzcyR9U9IZ4elXSzpN0nZJP5T0qhraDACYB3EbAGY3MUl295eNeeh5I57rks4t05AydTcDTEzfDCalB9LQRNyeJWYPdD12V7GOmkBsB0aL5op7OzacLql8DVcf6iZ3hnU0nzZr4PqwDQBkZo3ZA12OGzsij9lFdXkbAfMpW5MMAAAAdBZJMgAAAJBDkgwAAADkRFOTHLv56sYY1ABUgwueoCrEbKB+XY/ZHEmuAIMagGpwwRM0gf4CVKPrMZskGQAAAMghSQYAAAByOlWTPLd+cxI1Ll02ac5PU7G5QwH0U9kaR2BWXa+vxfSiO5K8bMmi0q+dtnPP8rdiFftnSuHqU11Vtm/E3qfQrqr7R98SZPaveMRYX1umf9CnqhPdkeQt553a2BWIJv3yS+FKSHn5z5TiZ0A9ONKBOjQZs/PyVyFNMd4RszEf4na7ojuSDAAAALSNJBkAAADIia7cogkMDEHXMOAEfURpAlJWJm4Ts5vVyyPJJMjomhgHnAAAxisTf4nZzeplkgwAAADMhyQZAAAAyCFJBgAAAHKiHLi3bMmi0nU3q9ZtOuiKbqvWbar9AhYxT9w9y7qsS36wTVevwle273V1fRRRtr/GvA/2QYxxZj4x95dY1+Vw3O5zjMLBuh6zo0ySByM3y4xcziclsybI+cnqU1NkFGzbI8S7ehW+sp+rq+ujCEZtp2mWmF01Ynb9+hyjcLCux2zKLQAAAIAckmQAAAAgJ8pyi1kNTkWlUvOCcqcPuzypOvV/6BNidnqI2QfLr48uf9Y+6fSR5BgHP6A6fdm+1P+hL/qyT/dVn7Zvnz5rl3U6SQYAAADKIEkGAAAAckiSAQAAgJyoB+61Pal6XwaRlF3PpvbrZScNHunL4ImuT+iONBCzm9H2ep4FMfuAMtuxL308FlEnycM7SlOTp6c+EX0ZVQWktie4HyXVL5Jp9eVLBXEjZjejyv09trjdl5gtEbdTQLkFAAAAkEOSDAAAAOREXW4xrKkarLn1mwv9nT7VTRWVcp1cHfq0PoruN3nsR91FzE5DDGNLYkHMnqxv+1EySfKW805tpHaqaKfpy440jVE7Tmz1bk3q0/oouz+wH3UXMTsN+at5djVGFUHMru91qaLcAgAAAMghSQYAAABySJIBAACAnGRqkqX6i+r7VLTflBjWaUyT15dZH1ZTW4DUxRBfUL3UY/bgdUjfxCTZzC6W9CJJd7v7E8Oy8yW9WtI94Wlvdverw2NvknSOpPslvc7dr6mqsVUMBJk08XxXi/Tbkg9kMa7fJr9k+zQqGO2JKW7PipiNPGI2mlKk3OISSS8Ysfyd7n5S+DcItCdKOkvSE8Jr/srMFlTVWABAIZeIuA0AM5mYJLv7ZyV9p+D7rZF0ubv/2N13SNou6ekztA8AMCXiNgDMbpaa5Nea2SskbZH0u+7+XUkrJF0/9Jw7w7JDmNlaSWsl6fjjj5+hGe0pcpqvbxNvAykrO8H+wqNXPaWG5tShdNwmZgOITd0xu+zsFhdK+klJJ0naLekd076Bu2909zl3nzvqqKNKNiN+DCoB0lF2f7XDFqQwCHqmuE3MBhCbumN2qSTZ3e9y9/vd/QFJ79GBU3O7JB039NRjwzIAQIuI2wAwnVJJspktH7r7Ekm3hNtXSjrLzBab2SpJqyV9abYmAgBmRdwGgOkUmQLu/ZJOkbTMzO6U9BZJp5jZSZJc0k5Jr5Ekd7/VzK6QdJuk/ZLOdff762k6AGAU4jYAzG5ikuzuLxux+KJ5nv82SW+bpVHz4YIiaYt1/Q4P6GHgzvSYcD8uscTtKvb3ufWb590fY40pqUplfQ7HbJO0Y8J82jgYMbuYFAabHGQQLOuaQL5ocsQE9uUUWb9tr9sUviBiw48KjFLFBaAm7Y/E7GqNWp+xrztvuwEJImYXU3Z2CwAAAKCzSJIBAACAnOTKLepQdDJqalXRB6vWbZr69CX7BppEzAYOKHNBDeq4i+FIsorXoFKrij4oU9/HvoEmEbOBA8r0c+q4iyFJBgAAAHJIkgEAAIAckmQAAAAgJ9mBe2Unwl65blMlgzmK/v2+TbxdhRgms580LygDgrqpbN/zB+7fX0NzOqGK/bmKuE3MLi+GmDzJpJjNQLVuqjtmJ5skz3JRkSp2dhKk+nDBEbSl7H5tF7zoxoqb0hlVXQBq1n2OmF1eCjF5EgaqdVPdMZtyCwAAACCHJBkAAADISbbcYlZV1Sajv6hbBpo1t34z+xRKK1ISQtzGsF4fSaauFHWifwHVYp9C3ehjGNbrJBkAAAAYhSQZAAAAyCFJBgAAAHKSH7hXxSTnTDKfnhQmt0+Vafo5Rdk3MEmV++w078Vgv2YQk9szy8XV8u/DvnKw5JPkLeedOvMk5nSK9BTdZrFPcB8jrkqFOlR1UZHBexV9HxK3ZnThgiOpyq/7suuZfeVQlFsAAAAAOSTJAAAAQE7y5RZSuRrKAWpy0Ka59ZtLneKinyJVs8RriVP2KZt128di1bpNpcZtELPT04kkeceG0ysLnNTkdMvOXH1tbF+wZfsb/RSpGtS8x7Yvon5Fxjuk0C/KJPrE7DRRbgEAAADkkCQDAAAAOSTJAAAAQE4napKBGFnbDQAg6UCdK4OngPmVqQnv8n7FkWSgJl0YxQ10CYOngOp1eb8iSQYAAABySJIBAACAnM7UJHdlknLUq+l+ksKcn0AbiNkYJYZ+kXrcjmEddkVnkmQmqEcRXZnMHkhdlReBQnfkYzR9ZHpFvudGYV0finILAAAAIIckGQAAAMghSQYAAAByoqhJvnnX3pG1MG1MUL1syaLCz51bv7nU/IBdnngbzZimnwJVGxezpfjjW9m4nRf75wQwu4lJspkdJ+kySccoGzC50d3/zMyOlPQBSSsl7ZR0hrt/18xM0p9JOk3SDyW90t2/UqZxTU5QvbNEoXvZ9nV54m1Mr0zfA8ZpM2ZL8ce3qtoX++cEMLsi5Rb7Jf2uu58o6RmSzjWzEyWtk3Stu6+WdG24L0kvlLQ6/Fsr6cLKWw0AGIeYDQAVmJgku/vuwVEFd79X0jZJKyStkXRpeNqlkl4cbq+RdJlnrpd0hJktr7zlAIBDELMBoBpT1SSb2UpJT5X0RUnHuPvu8NC3lJ3ak7JgfMfQy+4My3YPLZOZrVV21EKLHnPC2L9ZtH5s1vqwles2UWM2D+qv6zNpbkpT+Xkv0W9txGypurpflFfnNiCuoy8Kz25hZkskfVjS6939+8OPubtrygu8uPtGd59z97n5nld0Jx88b5YBTQT18fpUfx3boDiunIQy2orZUvP7fdF9NrZ9u051boMmt2+ftlnbyq7rLm+jQkeSzWyhsmD7Pnf/SFh8l5ktd/fd4dTc3WH5LknHDb382LCsEYNft1w5BmUVOUJC/0LMUorZUjZ4dZZ9apqjmuy7aanqiDXbfTLODhxq4pHkMPL5Iknb3P1Phh66UtLZ4fbZkj4+tPwVlnmGpL1Dp/gAADUiZgNANYocSX6WpJdLutnMbgjL3ixpg6QrzOwcSd+UdEZ47GplUwltVzad0KsqbTEAYD7EbACowMQk2d0/p2zs0CjPG/F8l3TujO2a2bIli0rVTeVPyTBAAaOU7V9NYrBlP6UYs+fWb57p9cNxm/7bDL4rp1OkbpeYHZ8orrhXhy3nnVpJDVLsiRDaMSogxVbz1qfBlkhblX2O/tsO1vuhpr1QFDE7PoVntwAAAAD6giQZAAAAyCFJBgAAAHKirkmeZXDUqgbqQ8u2r8sTb/fFqnWbGr3IR77emavwIVYxDGptYnwAcTweTcfjoob7YYoxu+x6TfGzjhNFkvykFUu1ZcwKLRvsmthhGE3aX20H5Lb/PvptvpgtxTeINW/aAVWIWwrxMIU25pVtc4qfdRzKLQAAAIAckmQAAAAgJ4pyi227vz/z6bmu1YfFNqk49dfxSb3eDemqImZ3TVdiNupT5z4z6b254Eg5USTJ+x8oXsHSl1qy2CYVZ+eKW5dqwBC/aWJ2X3QlZvPjp5v4wVQO5RYAAABADkkyAAAAkEOSDAAAAOREUZMcszoGnpUd4AHMZ9QFR8pUjjLYEjhYn2I2g7SBA0iSx6hzgGBfgi1GG+5bdQ6ScfVnoCviNam6KrW4AAAWtElEQVQPpjBQrE8xm0Hah8r34ab7bFPfGTgU5RYAAABADkkyAAAAkNPZcotZJ1KfW7+Z004Yq6q6vbon/M+fmmNCecSmjYteDO8X7BPpq7MPjaq1LjveAwfEdvGdcTqbJI9aidPU8vSpBg3Tq2onnbWfTot+jdjMty81UX/JPpG+pn/k5K9uSp3w9GK7+M44lFsAAAAAOSTJAAAAQA5JMgAAAJCTVE1y05OVDwbvFS0wj2EACBO6Y5LUBy2lMuADmUn1mn3fLk3GbPadNFVR89x2brBq3aYkBzsmkyS3cVGEQTApGlSaKCjn4hCoUoqDllIZ8IFiur5dYorZ7DvdFVM/GyXFBFmi3AIAAAA4BEkyAAAAkJNEuUWbtTR1zH9Y1YUogCpUccERah1RpbL9qS7E7P6KrS9Wpc7P1aV+H0WS/KQVS7Ul8nqaKpEUIGZlAie1jv1SJGbPcoAhtn5BzO6v2PpiVWb5XLHXP1eJcgsAAAAghyQZAAAAyCFJBgAAAHKiqEkeqPuiHaZ05+obhcFS3VR2kFCVUr/gCJrRhUFNKVzMo+x3F/tuf/Tloj1NDwqMKkmu+6IdO0KxeR0zVrSBwVLdVCaQ1dmn6S8Ypwt9o8nEoez6KntwpwvbJ2WTBrg1mYuk2hfaHiRIuQUAAACQQ5IMAAAA5ERVblG3uuvnujSBNtoTY51nFRccqQoXduiWNrdL3/oE+06zulLaOa2qvsNi6HcTk2QzO07SZZKOUVYatdHd/8zMzpf0akn3hKe+2d2vDq95k6RzJN0v6XXufk0NbZ9a1YlH27Uy6KbYEuRR2mxjFwaf1CmFmB1L7OxbX+rb50U7unShkiJHkvdL+l13/4qZPULSVjPbHB57p7u/ffjJZnaipLMkPUHST0j6lJk9zt3vr7LhAICRiNkAUIGJNcnuvtvdvxJu3ytpm6QV87xkjaTL3f3H7r5D0nZJT6+isQCA+RGzAaAaUw3cM7OVkp4q6Yth0WvN7CYzu9jMHhWWrZB0x9DL7tSIAG1ma81si5ltueeee/IPAwBmRMwGgPIKD9wzsyWSPizp9e7+fTO7UNIfKqt5+0NJ75D034u+n7tvlLRRkubm5pK7xkcMBeVAm6oYlMJ+VJ+YY3aMg1PRvBQuiBXDxZ2qQsyeXqEk2cwWKgu273P3j0iSu9819Ph7JF0V7u6SdNzQy48NyzojtsJyQIpr4voB9pV2xB6zu5J0YDYpXBCrSDI+S2yNMW4P63sMn1huYWYm6SJJ29z9T4aWLx962ksk3RJuXynpLDNbbGarJK2W9KXqmgwAGIeYDQDVKHIk+VmSXi7pZjO7ISx7s6SXmdlJyk7d7ZT0Gkly91vN7ApJtykbZX0uo6QBoDHEbACowMQk2d0/J8lGPHT1PK95m6S3TduYorU/42pimqpzG5z+aPOCCoO/z8TwsyvTb9re9qmY9lQh63V2McbsWV/TFWU/uyn7ZVPm73XVtLHFJO3oeelAGUXWc5fjdlRX3Jt1JTcdeNsO9F3tlE0rsx3b3vZdxXpNS9kY1HadZVuI2e1JbnaAhHQ5bk81BRwAAADQByTJAAAAQA5JMgAAAJATVU1yikbV1nW5iB31KzO4Z9QorSreN1UpXKQA7SN+owp1DqJPMW7PMkC1zIDMMvXmC49e9ZQiz+ttkjw8QXbVg0hS69CIS11f0EXftwuDqlK4SAGmN4jbdfZR+gCmVeePqrovZlKHJgf0lh2QaYctKJT/Um4BAAAA5JAkAwAAADm9Lbeo28p1m6htAxI06ZQf+/X0qrrQU2ynlfuqqpr/svWkXdXUBdFQHElyjejsQPewX0+PddYtVdX8l70CXld/LLGfxIdyCwAAACCHJBkAAADIIUkGAAAAcjpVkzzNBNYMrMNA2YnPGeBVD9YrgPnUFbMH7018mV5X43ankuTBBiha1E+RPKT6Jmunf9WD9Qr0W50X2CC+1CPV9Uq5BQAAAJBDkgwAAADkdKrcooyuzreINNQ5eXyqNWBV4CIFKMvabgCit2rdJuJLDRqN2wX/UO+TZKBNddZplX3vnSUm+I/tx2aRixTE1mbEoewFLtAfMSbIw3E71dhWZ9zOf6/ZBS/aWuR1lFsAAAAAOSTJAAAAQA5JMgAAAJBDTXKNli1Z1HYTgEaUndy/TWXbzH49vVT6B9u2GPaduHRlvdY1kH2W9UOSXIEyA52ALhk1i0bsg0f6OvNHGyat6xj6CnG8OPad+vS5HzY52Lwoyi0AAACAHJJkAAAAIIdyiwqsXLep1xduAJo26fS8iflu21K0rpCYWa2y9ZxsBzQlhrKqaXXySHIbRewpDEpBeWX6VFcGU5TV5uePcbL/vigaC4ef1/d9pQplv4O6+t1Vtk91oS924TPEopNHklMcRIS4caRlevl1xj6IcYjZqFqfY3bZz84+d6hOHkkGAAAAZkGSDAAAAOR0stxilFQms0e/1Nkvi9SlMdgHsWo6Zs+t30yfRiGmesY9ELPLyZeJVPlZe5Mkj1th1OCgTW0HLQb7IFZNx2z6NIpqc+YcYvZkVX5Wyi0AAACAHJJkAAAAIIckGQAAAMiZWJNsZg+R9FlJi8PzP+TubzGzVZIul/RoSVslvdzd95nZYkmXSTpZ0rclnenuO2tqf1RG1cpx5S8g0/RArPz+WGZfXLVuU6EBOosec8LJU71xjboYs4tuh1kN+gxxG01oql8X1aUJDiaNXSgas4sM3PuxpOe6+31mtlDS58zsHyW9QdI73f1yM3u3pHMkXRj+/667n2BmZ0m6QNKZRRrTRTHtAECbigxSrHMgbZl9MdH9t3Mxu+ntkOh2R2Ji62dlBpJ3ffKDieUWnrkv3F0Y/rmk50r6UFh+qaQXh9trwn2Fx59nZlZZiwEAYxGzAaAahWqSzWyBmd0g6W5JmyX9m6Tvufv+8JQ7Ja0It1dIukOSwuN7lZ3eAwA0gJgNALMrNE+yu98v6SQzO0LSRyU9ftY/bGZrJa2VpOOPP37Wtyus7ETcQEyYUL6c4VODXa47bTNmF61rLHLhBBRXtp6U7dAMco80TXUxEXf/npldJ+mZko4ws8PDkYdjJe0KT9sl6ThJd5rZ4ZKWKhsMkn+vjZI2StLc3FxjpTl0UnQBE8rPLrZ6wDq0EbP7/COsTaz3uBF70zSx3MLMjgpHI2RmD5V0qqRtkq6T9NLwtLMlfTzcvjLcV3j80+7eh+8jAGgdMRsAqlHkSPJySZea2QJlSfUV7n6Vmd0m6XIzWy/pq5IuCs+/SNJ7zWy7pO9IOquGdgMARiNmA0AFJibJ7n6TpKeOWP4NSU8fsfxHkn65ktYBAKZCzAaAakxVkwygW9oe7NP2YJb8HJ99H9gIoJtM089pPGpwc9sxu2kkyUCPtZ0QxhZsY2sPAOTtLDErT5mLfowamNC3GFlonmQAAACgT0iSAQAAgBzKLWrGRO1AWsqclkS3cE1uxIqcolkkyTUqUzcEID6T9mW74EVbG2oKGsAk0YgJuUR7KLcAAAAAckiSAQAAgBySZAAAACCndzXJZS+eUMbKdZtGTsYNzKrti4D0zaTBfIsec8LJDTWll5qM20AdUo3ZfR/I3LskuejFE6rqGAwAQR3avggI0KRx/b3vX+BIBzE7TZRbAAAAADkkyQAAAEBO78ot5jO3fjN1b0ABZfaVZUsWHXLKkVrTfinab0b1lVFMlLTFquz3adFtj+mwr5RDkjyEL2ugmDL7yqjXFPkypO60O4r2m6LP40s/XmW/T/kerkeV+8qki5t0KWZTbgEAAADkkCQDAAAAOSTJAAAAQA5JMoCotT2ZPgCguC7FbAbuAYhaVSPduzSYBABiNSpmpxp/OZIMAAAA5JAkAwAAADkkyQAAAJioS/XGRVCTDAAAgENMunBI13EkGQAAAMghSQYAAABySJIBAACAHHP3ttsgM7tH0g8k7WmzHQuPXvUUO2xB5XXa+761fWvV7zmDZWp5PZeUYrs72+Yy+4o/cP/+/7h7x42lWzZeoTYveswJJ9fwtyVJ+/ferft/uNfqev/YlI3ZRftN0b5SV8yWoorbKcYRLTx61Ul22IIF076uxjhRRIrruraYLZf23VXNflBn/C2jaMyOIkmWJDPb4u5zbbdjGrS5OSm2mzY3I8U2d0GK6502NyfFdtPm5qTSbsotAAAAgBySZAAAACAnpiR5Y9sNKIE2NyfFdtPmZqTY5i5Icb3T5uak2G7a3Jwk2h1NTTIwzMx+UdKJ7r6hxGt3Sppz9z1mdp+7L6m8gQDQMjP7jKQ3uvuWGd/nJEk/4e5Xj3l8TtIr3P11M/4dk3StpBe7+/dLvP71kja6+w9L/v1TJO1z938J939d0g/d/bJ5XnO+pPvc/e1mdomkq9z9Q2X+fszMbJGkT0l6rrvvb7s9sYjpSDIgSTKzw939yjIJcteYGVfFBFpkmcq/KyPbt0+SdNqoB0I83jJrghycJunGMgly8HpJD5vh758i6WcGd9z93fMlyF0zX59z933KfsCc2VyL4keSjEOY2Uoz22Zm7zGzW83sk2b20PDYZ8JRBZnZsnDUVmb2SjP7mJltNrOdZvZaM3uDmX3VzK43syPD837SzD5hZlvN7J/N7PFh+SVm9m4z+6KkPwrv9xfhsWPM7KNmdmP49zNh+cfC+9xqZmsnfKbLzOzFQ/ffZ2Zrcs85JXy+D5nZv4bnWHjseeGz3GxmF5vZ4rB8p5n9gZl9JTz2+BF/e87Mbgj/bjYzn3JdHBk+601hXT556o0KoLAQA283s8sk3SLpODN7vpl9IezrHzSzJeG5G8zstrB/vj0sO8rMPmxmXw7/nhWWn29m7zWzz0t6b9ifnzD0dz8T4sXDQ5z5Uog7a8LjDzWzy0N8/qikh45p/8lm9k8htlxjZsuH3v+C8L5fM7PnWHYE8a2Szgwx6swR7TzFzK4K7zGubU8Iy24I62L1iKb9qqSPD63jQZzdFuLuw8Jjh8RbM3udpJ+QdJ2ZXReeN26bHBKXzWylpF+X9Duhjc8Jn/ON4TWvDtvqxrDtxibjZvZcM/vY0P1Tw/bIP2/k98O4mB7ac3HYTt8In3nU37/aDnyn7DWzs81sgZn9cfgMN5nZa8JzT7Hs++VKSbeFZW8ws1vCv9cPvfXHwjbCgLvzj38H/ZO0UtJ+SSeF+1dI+rVw+zPKShmkbH7GneH2KyVtl/QISUdJ2ivp18Nj75T0+nD7Wkmrw+2flvTpcPsSSVdJWjD0fn8Rbn9g6PULJC0Nt48M/z9U2RfZo8P9nZKWhdv3hf//i6SPhdtLJe2QdHjuc58S2n2ssh+QX5D0bEkPkXSHpMeF51021J6dkn4r3P5NSX8zYd3+saQ/nnJd/Lmkt4Tbz5V0Q9t9hH/86/K/EAMfkPSMcH+ZpM9Keni4/3uS/rekR0u6XQdKF48I//+9pGeH28dL2hZuny9pq6SHhvu/I+kPwu3lkm4Pt/+PDsTcIyR9TdLDJb1B0sVh+ZOVxem5XNsXSvoXSUeF+2cOveYzkt4Rbp8m6VPh9isV4u2Ydp6irMxgvrb9uaRfDcsXDV6ba9s3JT1iaB27pGeF+xdLeqMmx9tl822ToecdEpfD53pj7nO+Mdx+9NDy9UOvH37OJZJeKskk/evQOv57Sf91xOcd146RMT38rX+RtDh8vm9LWjhPPz1Z0k3KvtPWSjovLF8saYukVWHb/UDSqqHX3By22RJJt0p6anhsgaR72t7/YvoX0+kexGWHu98Qbm9VFtAmuc7d75V0r5ntlfQPYfnNkp4cfuX/jKQPmj04h/fiodd/0N3vH/G+z5X0CkkKj+8Ny19nZi8Jt4+TtFpZUDmEu/+Tmf2VmR0l6ZckfdhH1119yd3vlCQzu0HZ575X2fr4WnjOpZLOlfSn4f5Hwv9bJf23UX8/vN+Zkp4m6flTrotnhzbL3T9tZo82s0d6+VOWACb7prtfH24/Q9KJkj4f9tdFyn5E75X0I0kXhSOtV4Xn/7ykE4f27UcOjnJKutLd/1+4fYWkT0p6i6QzJA1qXZ8v6RcHRzmVJY7HS/pZSe+SJHe/ycxuGtHu/yzpiZI2h7+/QNLuoceH49XKeT7/cDuHjWvbFyT9vpkdK+kj7v71Ea89MnxHDNzh7p8Pt/9O0uskbdb88XZg3DYZ9TnHxuUhTzSz9coS/yWSrhn3RHd3M3uvpF8zs7+V9EyF76gRRrVjZEwPj21y9x9L+rGZ3S3pGEl35t/UzJZJeq+kM9x9r5k9X9n37EvDU5Yq+07cp+x7bcfQ3/6ou/8gvM9HJD1H0lfd/X4z22dmj8htp94iScY4Px66fb8OnNbbrwNlOg+Z5zUPDN1/QFlfO0zS99z9pDF/8wdFG2fZAIyfl/RMd/+hZQNY8u3Ju0zSr0k6S9Krxjwn/7mL7COD1zz4fDO7Rllw2+Lu/8PMnqjsKMHPhkBU2boAUIvhfdAkbXb3l+WfZGZPl/Q8ZUcYX6vsR/1hyo5C/yj33IPe1913mdm3w+n2M5WVAwz+3i+5++0jXj+JSbrV3Z855vFD4tUY42LQyLZJ2mZZidjpkq42s9e4+6dzz9lvZoe5+wPhfn7mgGlmEhi7TYKin3PgEmUDCm80s1cqOwI7n79VdiDoR8oOaowb7DZtOw75DjKzcyW9Oiw7TdJdki6X9FZ3vyUsN2VHrQ9K7sN35TTfJ4uVfSaImmRMb6ey0zVS9qVQWDjyucPMfll6cEDMUwq89FpJvxFes8DMlir7lfzdkCA/XtlRhUkuUTbwQ+5+2xRNv13SSjM7Idx/uaR/mu8F7v4L7n5SSJCPkPR+ZaPD7wmPT7Mu/lmhTiwEvD0cRQYadb2kZw1igGV1uY8LR4eXejYrxO9IGuzDn5T0W4MXWzZ7xDgfkPQ/w/sMjgxfI+m3zB4cE/HUsPyzkn4lLHuispKLvNslHWVmzwzPW2hDdc9j3KusVK6IkW0zs8dK+oa7v0tZ3fG4tj126P7xg3Yq+1yf0/zxdridI7fJhLbP9zkfIWm3mS1Ugbpcd/93Sf8u6TxlCfM0porp7v6X4fvkpPB3N0i6yd0vH3raNZJ+I7RfoX8+fMzffrGZPSw8/pKwTGb26NCW/5jy83QWSTKm9XZlO+JXldVMTetXJZ1jZjcqq4VaM+H5kvTbkn7OzG5WdsrqREmfUPYLe5uygHH9PK+XJLn7XZK2acqAFo4GvUpZacTNyo6Mv3uKt1gj6T9Jes9gsEVYXnRdnC/p5HBqdYOks6dpP4DZhB+3r5T0/rAffkHS45UlVleFZZ9TVjMsZWUDc2EA1W06cIR4lA8pO7t1xdCyP1RWW3yTmd0a7kvShZKWhLj3VmXxMN/WfcoOYFwQYssNGprRYYzrlJWH3BDKwuYzrm1nSLolxLcnKjtzl7dJBx+hvV3SueHzPErShRPi7UZJnzCz6+bZJvP5B0kvCZ/zObnH/pekL0r6vLJ64yLep6xkZFvB5w+cr9li+huVle0NBu/9oqS/UTYw7ytmdoukv9aII9fu/hVlB4y+pOzz/o27fzU8/HPKthEC5klGb1g2WvlmSU9z972Tng8AqI5ls2xc5u6nWjbbxFXu/sR2W1WeZTMwfdXdL2q7LVUI9cnrhurBe48jyegFM/t5ZUeR/5wEGQCa5+67lZ1Re+TEJ0fOzLYqKyn5u7bbUgXLpgL8GAnywTiSDAAAAORwJBkAAADIIUkGAAAAckiSAQAAgBySZAAAACCHJBkAAADI+f+jkFK+FgHz8gAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 864x864 with 2 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "try:\n",
    "    import scipy.sparse as sp\n",
    "    import matplotlib.pylab as plt\n",
    "    plt.rcParams['figure.figsize'] = (12, 12)\n",
    "    rows,cols,vals = a.mat.COO()\n",
    "    A1 = sp.csr_matrix((vals,(rows,cols)))\n",
    "    A2 = A1.copy()\n",
    "    A2.data[:] = 1\n",
    "    fig = plt.figure()\n",
    "    ax1 = fig.add_subplot(121); ax2 = fig.add_subplot(122)\n",
    "    ax1.set_xlabel(\"numerically non-zero\")\n",
    "    ax1.spy(A1)\n",
    "    ax2.set_xlabel(\"reserved entries (potentially non-zero)\")\n",
    "    ax2.spy(A2)\n",
    "    plt.show()\n",
    "except ImportError:\n",
    "    pass"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "#### Remark 2: Dof numbering of higher order FESpaces "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "-"
    }
   },
   "source": [
    "* In `NGSolve` `FESpace`s typically have a numbering where the first block of dofs corresponds to a low order subspace (which is convenient for iterative solvers). \n",
    "* For L2 this means that the first dofs correspond to the constants on elements. \n",
    "\n",
    "* You can turn this behavior off for some spaces, e.g. for L2 by adding the flag `all_dofs_together`.\n",
    "\n",
    "We demonstrate this in the next comparison:\n",
    "\n",
    "\n",
    "(you need matplotlib/scipy for this block - but can ignore it otherwise)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAA2oAAAEkCAYAAABe5sI+AAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMi4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvhp/UCwAAIABJREFUeJzsvX1wHdd55vm8Er/Fb8giKIFYwKGYyA5iioQ+qU0cK2ZMxiONKrLKTkprU06xdsvJOg4zsZzsjsytfNgzZhRtZcoe1Vj+mHVM24w1YsXihhytPVnRS3IEkg4CMqYYgiZBEZQFUAQsgiQgnv3j3r44t3H63tN9++N09/OrQqFvo+/t093nHvTb53mfV5RSIIQQQgghhBDiDjdk3QBCCCGEEEIIIfUwUCOEEEIIIYQQx2CgRgghhBBCCCGOwUCNEEIIIYQQQhyDgRohhBBCCCGEOAYDNUIIIYQQQghxjNQDNRH5gIj8WEROisiTae+/GSLynIi8LiL/pK1bLiL7ROTV6u9lWbbRQ0RWicj3ReSYiAyIyCer611t7zwROSQiP6q2d3t1fbeIHKz2iW+JyJys26ojIjeKyBER+bvqa2fbKyKnRaRfRI6KyCvVdU72B1dIa0zK8vuadh8WkaUisktE/llEjovIfSkd56eq5/afROSb1TEn1mMN8z9CKvyf1X3/o4isi3Gf/756fv9RRJ4XkaXa3z5T3eePReTX4zxW7W/bRESJyM3V14kda3X971WPd0BE/p22PpZjdZU0xqcyjU3VfaQ+PqUxNlX3U4rxqVRjk1IqtR8ANwL4FwDvBDAHwI8AvCvNNli08ZcBrAPwT9q6fwfgyerykwA+n3U7q21ZCWBddXkRgBMA3uVwewXAwurybAAHAdwL4NsAPlxd/yUA/0vWbfW1+w8A/A2Av6u+dra9AE4DuNm3zsn+4MJPmmNSlt/XtPswgK8B+J3q8hwAS5M+TgC3ARgEMF87xo/Ffaxh/kcA2AxgT3XsuxfAwRj3uRHArOry57V9vqvaj+cC6K727xvj2m91/SoAfw/gJ954k/Cx/iqA/wpgbvX1LXEfq4s/aY1PZRqbqp+b6viU1thU/ZxSjE9lGpti7fwWB3kfgL/XXn8GwGfSbINlO7t8F+LHAFZWl1cC+HHWbQxo9wsA3p+H9gJYAOAwgHsAvKF9oev6SNY/ADoAvATgfQD+rvpFd7m9pzEzUHO+P2R4vjIbk9L6vqbdhwEsQeWmRHzrkz7O2wCcBbAcwKzqsf56Esdq+z8CwH8E8BHTdq3u0/e3RwB8o7pc14dRuWm5L65jra7bBeA9+niT5LGickP7a4btYj1W136yGp+KOjZVPzP18SnNsan6WaUYn8oyNqUtffQ6q8dQdZ3rrFBKna8uDwNYkWVjTIhIF4A7UZmlcra9VZnDUQCvA9iHylOGN5VSU9VNXOsTfwXgjwBcr75ug9vtVQD2ikifiGytrnO2PzhAJmNSyt/XtPtwN4CfAvhKVdL0n0TkJiR8nEqpcwC+AOAMgPMALgHoQzrf16BjS6t/PYHKE+PE9ykiDwM4p5T6ke9PSe53DYD/sSoT+28iclcK+3SB1I+v4GMTkMH4lPHYBJRkfCrq2EQzkZCoSmissm6HjogsBPC3AH5fKTWm/8219iql3lZKrUXlSdrdAH4h4yYFIiIfBPC6Uqov67aE4AGl1DoAmwB8QkR+Wf+ja/2hjKT5fc2oD89CRR7yRaXUnQDeQkVuUyOJfljNu3gYlRuxWwHcBOADce7DhrS/YyLyJwCmAHwjhX0tAPDHAP5t0vvyMQuV2Yh7AfwbAN8WEUm5DYWnBGMTkMH45MrYBBR3fCry2JR2oHYOFf2oR0d1netcEJGVAFD9/XrG7akhIrNRGVi/oZT6bnW1s+31UEq9CeD7qEz3LxWRWdU/udQnNgB4SEROA9iJijzjGbjbXu/JHZRSrwN4HpVg2Pn+kCGpjkkZfF+z6MNDAIaUUgerr3ehcmOUdD/8NQCDSqmfKqUmAXwXleNP4/sadGyJ9i8R+RiADwL47eoNWNL7/DlUbjZ/VO1THQAOi0h7wvsdAvBdVeEQKjMwNye8TxdI7fhKMjYB2YxPWY5NQDnGp8KOTWkHav8dwO1ScbqZA+DDAHan3IYo7Abw0eryR1HRb2dONWr/MoDjSqm/1P7kanvf4Tn/iMh8VDTwx1EJ2B6tbuZMe5VSn1FKdSilulDpq/+PUuq34Wh7ReQmEVnkLaOSzPtPcLQ/OEJqY1IW39cs+rBSahjAWRH5+eqqBwEcQ/L98AyAe0VkQfVce/tN4/sadGy7AfxPVdexewFc0iRILSEiH0BFNvaQUuqyry0fFpG5ItIN4HYAh+LYp1KqXyl1i1Kqq9qnhlAxoRhGgscK4L+gkrQPEVmDigHEG0jwWB0hlfGpLGNTdb9ZjE9Zjk1ACcanQo9NURLbWvlBxX3lBCq5SX+S9v4t2vdNVDTEk6hc6I+jopt+CcCrqLi7LM+6ndW2PoDKFPY/Ajha/dnscHt/CcCRanv/CcC/ra5/Z7UDnwTwHVTdc1z6AfBeTLtSOdneart+VP0Z8L5frvYHV37SGpOy/r6m2YcBrAXwSvVY/wuAZWkcJ4DtAP65Or78Z1Qct2I91jD/I1AxR/gP1b7VD6A3xn2eRCUHwutLX9K2/5PqPn8MYFOcx+r7+2lMJ+wneaxzAPxf1et6GMD74j5WV3/SGJ/KNDZV95H6+JTG2FTdTynGpzKNTVL9MEIIIYQQQgghjkAzEUIIIYQQQghxDAZqhBBCCCGEEOIYDNQIIYQQQgghxDEYqBFCCCGEEEKIYzBQI4QQQgghhBDHyCRQE5GtWew3KmxvsrC9yZK39saFiHxARH4sIidF5EmL7VM/T1ldm7IcK88v9+kieRibstpvWfaZ1X7Lss+s9pvEPlsK1MIONhp5G2jZ3mRhe5Mlb+1tGRG5EZW6KZsAvAvAR0TkXU3elsV5yuralOVYeX65T6fI0diU1X7Lss+s9luWfWa1X3cCtYiDDSGEpMHdAE4qpU4ppa4B2Ang4YzbRAghHJsIIdZELngtIvcB+KxS6terrz8DAEqpvwh6z80336y6urrw05/+FO94xzsi7TcL2N5kYXuTJen29vX1vaGUcuqEiMijAD6glPqd6uvHAdyjlPpdbZutqD79mrtg4Xq1uB23LJqLFYvnpdbOsNfmwtgVvD5+teV2JtUnvPYBmNHGLL43WX1X49xvo3Oa1D5tcX2feR2bqutr49MN8xevX/KOW/HWtSkAjftBnLh+ffO8T9v96t//m+bMwlvXplq6/jy/7uzTdnya1UJ7bgNwVns9BOCeRm/o6uoC/vVfQH52DW8AuHnhHLzyv72/hSaQvNH7p/vwxs+uAeD1Lwoi8pOs2xAFpdSzAJ4FgF9au0795v/xnwEATzzQDQDY1TeER9d3YNXyBZm10c+hwRHs2HsC2zauwd3dbVk3ZwZnRy/juZcHAQCbetqxp38YQOWceudR30ZfT8wcGhzBUy8MAAC2P/xuJ6+7q+R1bALqx6fl/8Md6ulvfg87Dw0BYD8oE/r3/8N3d+DF/gvOjv8kHLbjUyuBmm1Dak+FOjs7IdWbdAC1G3ZSHt7g9SfpcA7AKu11R3WdkTmzbsDi+bPxzEuvAgCOnLmIo0OXMDYxiSce6HYmaNt/cgQHB0ex/+SIk/+oVy1fgKceejcA4Ol9J/CVH56esc34lUnsOjx9KbztGcCZ2X9yBMeHx2vLLl53EopQYxMAvHVtCjsPDdX6wZ7+YfaDkrCnf7h23b0+wOtfLloJ1KwGG/2pUG9vr3qjhR0SQogl/x3A7SLSjcq49GEAv9XoDY+u7wAAjE1M4ujQpdr6XX1DtQDuU+9fk1Bz7fDa6P12mUfXd2BsYrL22gva1nYsMW6/q2+ots3i+bMzP9euoJ/HPFx30pTQYxMAQBJuFXEf9oFS0kqgFmmwuXnhnDrpGykXvP4kDZRSUyLyuwD+HsCNAJ5TSg3YvHd911IcOLUIQEW+t3LJfABu3CSvWr6gYQBzdvSyM7N/Ovo5fWTdrbg6dR1A5fx6MCAxo89ShoWzlO4RZWy6ZdFcbH/o3TUp8aaedmzfPVBb3n9yZMZ33rv241cmsWjebF7/nLKppx0HTo0AAD58V0X6qI+bpPhEDtSi3gj5c5JscpaY10QICYtS6kUAL9pu782c3dO9vE5ipONqIOTh0uyfPkN27Pz0OX2x/0Ld+dVvPk0w2IgOZyndJOzYtGLxPNzd3VaTu23fPVC7rgdOVaSxYxOTdQH9cy8PzpAeRw34SXZQ+khaylELO9iYsMlZYl5TceC1JK7izeJsWN1WF6DpN7oAnAmETLgkjdRnyHRjEX0ZgBbMjeHg4CiA+qCCwUZ0OEtZAiiHKw+81qUkcTMRQgjJA56s8OzoZSyePxuPru/A+UsTNdnJhtVtTskgdfSZPlcCGV2yd3b0cm39yiXza+sPDY7Uzu/mnhW1oGLD6umnxQw2wuGfgeQsSvGwkcN521y+NoUFc2ZRLpdTKH0kmQdqNjlLzGsqDryWxHV0+SAAowzy/KUJpySQLkkeTQTNiumOhs0kkeNXJvHcy4OBOTmkAmcgi4//e3NwcLRWtsP73oxfmaxtA8x0ivQe7ugKAkqL3cN0ren+Wi4yD9Rs8s2Yk1YceC2J6/jlgybnQl2m58KNsEuSRxNBs2JB8kjAfK79r104967BGcji4//e7Nh7ohased+NIHdVD+/hzoFTI0bJMXED/7X2HlKR8pB5oGbCbx4CgIYjhJBU0J0Vg2SQj9/XiXvf2ebMP0zX3SBtnAubSSIn376OtauW4rG7VvFmpQFRXSL1GRbOWLqN/xpv27gGO/aeqJMPe+6qQdLHDavbcODUSKDkmLiB/1pzJq18OBmoNTKcoOEIISQtdEnh2MS0lKjv9Js1c5GsgyAbXJVGBsn0dKez5w+/VltuWzh3xmcwwIgH0wyLS32FBLOnfxgHB0frxii9QLa3jX6Tb3oP3QQJcQ8nAzVCCHEBXVLoGTR4kiI9j83FIEjHVWmkjUxv9YqFuLNzGY6cuThD3vWp969hgBETuuspZyxzigQs276HEOIcTgZqJsMJGo4QQtJGlxQ+8UB3TQbpMTYxifVdS3FP93JsWN3m7OxaM2lkVgTJ9J54oLtuedXyBTPMD8YmJnFocATnLl7G2o4l2NyzAsC0fIszbeHQ+whnVdxGd/bc1NOO8SuTWNuxBI+suxU7Dw0BqDgE7jw0FCh99MYzv5kIKQf+PsQx0l0yD9Tiyi1jTlp5YX4iSQN/7tqRMxdxdOgSDpxaVCtCunj+bCdn11wNIINolmc1fe6nHdGuTl2vKwbLmbZwMLDND3oxa9N3AGgufWRgXm5MfchfNJ24QeaBmim3LEqOGikvzE8kabOrbwhHhy5VXmjSIVclhq7mqIVBD7xM594v4aKULxwMbHNK0HeAkkZiC/uK02QeqBFCSN4IspU/f2mi9tulGSxXA8gw6IGXXmPNtHxocKTO6l+HM0dmGNjmB38RZL/cUV9mwWtigoW080PmgZoptyxKjhopL8xPJGmzavkCPPFAN3b1DWHlkvkAKnW/PAnJf/vx67VZHxdmJaLkqLkml9SPYeWS+djVV7kh9fIGVy1fgD39w3XXAZgp6+HMkZk0pXB6fgyLLIcnqFB80DIws4i8adl/LVgUu7iwkHZ+yDxQM+UT2eQYsdYa8eB1JFmgywlrVCUk41ensLZjSWJ1idIIolyWSzYNtgwSsCNnLjY0H8kDRQlwgsoyEDuCZvRNy+NXJrFo3myMX5nErsPnAMD4IMNDz1Hy8pj2DgwHbkPyCQtp54fMA7WosNYaISRL/HJCz0FNt49P6gllGkGUy3LJIJmeycluU0977Zro12bu7Btz9xS5KAGOTVkGEkyjIshBy9t3D0x/QNhcNuYwFQ4W0s4PuQ3UCCEkS/wukEBFkrdt4xo89ULlpiip2RqbIKrVWTdXLf0BO5meLonctnENduw9gc09K2oBQh5n1OIKcLKemWvm6knixyavzZ+jtKmnHcfOj2Fzz4rAbQghyZLbQI211gghruCXQXoyoT39w5kZVrgsXUyKIEnk/pPTr/15GX/+veO4s3NZrqWEYSnKzBypYJNLtqd/2Gjdry/7LfxN35s8zUCTCsw1zDe5DdSi5CUxl4kQkgT+GS5v1gNAIsGSTRDmsnQxKYIkkUGOkZ4U8ujQpVwELHEFWJQeFgvTA4qG/cNS+mj63rC/5I/Q/YM4hZOBmo3pRxQzEVIMaApDXMMvE/RkXZ4kcmxisrYchwmId6Ptfa7ps1yWLiZFkCQyaP0XPvSemgRww+q2Wh6P/rQ5a5mgjh5gBbXXBkoPi8WG1W04cGqkTtrrl/RGkT6yKHYxsOkfxF2cDNRsTD9YFLu80BSG5IVVyxdg8fzZtdmvI2cu4ujQJYxNTNbs/aMEbfrn8slodPSA5el9J2qzVTq6Wx5QH4SnHcAFtZd9oNzs6R/GwcFRjE1MBsoYo0gfSTGw6R/EXZwM1AghpCjockivthrQeg5ZGaWNSaLPVgGoBUFrO5YYt886z4vyRTIDW3fGsK6PpBjwWucSJwM1G9MPFsUuLzSFIQAgIqsAfB3ACgAKwLNKqWdEZDmAbwHoAnAawGNKqYtZtdNjfddSHDi1CEBFYuQVyo56kx2XtNG1wtYuoF+rR9bdiqtT1wHUS8OyDpRcki/qZgVZmee4Rprjk6ksxRMPdM/YxqNR7TVPTmtbny1rOTBpjk3/IO7iZKBmk3PUbBtTHhPz2ooBrxOpMgVgm1LqsIgsAtAnIvsAfAzAS0qpz4nIkwCeBPDprBrpzZzd0728Tnaik2WwVEZ3SBP6DNmx88vr3CH166bfsJpwKactLZoWIC8nqY1PNrlkNrXXtu8eqH0Hgopi+wtkAyyA7TrMNcw3TgZqcWDKY2JeGyHFQSl1HsD56vK4iBwHcBuAhwG8t7rZ1wD8ABkGaibnNAB1sjkgGXfIMO0ru3xOnyELmk0AoAVzY0YHtawlkVkQ5LZZZvIyPgUSJI+kfI6QVClsoEYIKQ8i0gXgTgAHAayo3iQBwDAq0qPM8J5mnh29XCvCfP7SRM2BbcPqtpZlkHG0LypFkU7qMw6eQydQKWLurT80OFK7bkEOallLIrMg6yf2rksvXR6fdGycIb3ly9emsGDOLBbAJiRhChuomfKYmNdGSPEQkYUA/hbA7yulxkSmH/kqpZSIqID3bQWwFQA6OzsTb2ejotge5y9N5C7oKaJ0MmhWbP/JEWtJ5PiVSTz38iA29bQ7GTwUCZell66MT7okN6hP2jhDPn/4tTrpIwtglxcW0k6HwgZqpjwm5jYRUixEZDYqN0HfUEp9t7r6goisVEqdF5GVAF43vVcp9SyAZwGgt7fXeLMUJ42KYpvkdK7caDajiNLJoFmxIHkkYL6G/td5uaZ5w1XppUvjU30OpkWfDJA7rl6xEHd2LquZibhyrkn6sJB2OuQ2UGtW9NjGTISBG9Fh/8gXUnk0/WUAx5VSf6n9aTeAjwL4XPX3Cxk0bwa6PCxIBvn4fZ24951tid78xC1VLGJhbRtHxWaSyMm3r2PtqqV47K5VTgUPYXFdVghkL7004dr4tGF1G/YOVB4sbO5ZUVunYyN9fKx3lTPnmGQLC2mnQ24DtWZFj8OaiRDC/pE7NgB4HEC/iBytrvtjVG6Avi0iHwfwEwCPZdS+QHS5oF6EtO/0mzVzkaRyv4ooVUyLIEmkLhnTpWFtC+fO+Iw8BD46LssKHcep8ckka2RRbNIKLKSdDrkN1Agh5UYp9TKCPcgeTLMtYdHlgl7eiFdYWc9jSyKgKqJUMS1sjEI8adiRMxdxcHAUO/aeqAtw8hb4uCordB2nxycb50Y6PRJb2D8SJbeBWrOix7ZmIoR4sH+QtNClWl4xUv0GeGxiEuu7luKe7uXYsLot1tm1IkoV0yJIEqkXj/WS6f2J9mMTkzg0OIJzFy9jbceSGfIzV2faXJQVkvD4C16bgm6/9PH5w69h9S0L8dhdq2YUxQ4qeE2DifLAQtrpkNtALUr+kE0emx/mLZUHXluSBf7ctSNnLuLo0CUcOLWoJk9aPH92rLNrRbHUd4VmOW3T13TaNfLq1PU6+VkeZtpcDSZJcxoVvPYwSR/v7FyGu7vbjEWxPfTPfe7lQXzlh6exd2C4zh2SRbGLBx/ipENuA7Vm2OQbxbUNIYTEwa6+IRwdulR5oclJ4pYrMk8tefTAy3RN/XKhPEgM8xBMkphoVc5GORwhsVDYQI0QQvJGkP37+UsTge+JMjvGPLXk0QMvXTZmWj40OFJn9R8WvUZWkjKzPASTJDp+6eOL/RdmFLT2tgkqeL2ppx3Hzo9hc8+KmlMki2ITEp3CBmo2+UZxbUMIIXGwavkCPPFAN3b1DWHlkvkAKvW5PMnc2MTkDAlRlNkxf54apZDxo5/jlUvmY1df5abVy0dctXwB9vQP111fAA2vdRBBbpRxYyN1SitoJPHjL+h+cHB0RkFrfRtgZrH3P//ecRwdutTQCZB5bITYU9hAzSbfyCYnzQbmsRFC4kIPvGpUZUQHTo1g++6Bupsar5ZNK/VrKIVMlqaSQYMk8siZi7WZtiDjBg+9RlbWdYzSChpJ/Phn9E2zpt42Xp8cvzKJXYfPAUDdA4e3rk3V3jN+ZbLuM5jHRog9hQ3UotIoJ415bISQpPHLEj1XLc/m/fjweN0N8P6TI8Yn363sk8RLkGTQ5Jq2qae9dq11a38gOPDRZzla6QdxYFPCgLiJjeGIf5vtuwem/6g9cLhp7vTt5aJ5s807ZB4bIU1pGqiJyHMAPgjgdaXUL1bXLQfwLQBdAE4DeEwpdTG5ZhJCSHJcm7peu+HIWoLjlyV6y9s2rsFTL1TaqM+axBFkxWXZr8ve9BysrM9p1gRJBoPWb9u4Bjv2nsDmnhUYm5is5QMFzZa5FBw1c8AkxcKf1+blpenLzGMjJDo2M2pfBfDXAL6urXsSwEtKqc+JyJPV15+Ov3npY8pJYx4bIcXm4uVrscm19HwvALHlfumzJnpeyBMPdDsjL9Nlb8fOj9Vmg7KWwAXZyruaT+XNkgJomA+k9wEGRyQNvO+MLn30W/r7l/05aqb+nfVMMAkPcw3ToWmgppT6BxHp8q1+GMB7q8tfA/ADFCRQi6M+GyEkXyxbMAe/eX8XgNZnJPw5ZnHlfumzJgCczAMKcq3MepYnKEfM1Xwqv2Okd1MMwMlAmJQH/TsDAGs7lkz/sUH5CR2TI2rWYwQJj2lc5VgUP1Fz1FYopc5Xl4cBrIipPaliMgExmYnQKKQc0BSmvMyZdUNsMxJ+KeLYxCTGJiZxdvRyS08adUnZocGRmtyokXlE2m6ONjkuWRCUI6YHlhtWtzkpf9XPoX7dPVkkEN1AhAWsSVg84xpPjvvIultxdeo6AHvpIwslFwPPyCqOsYgE07KZiFJKiYgK+ruIbAWwFQA6Oztb3V2smExAopiJkGJAUxgSB/58r8XzZ+OZl16N9KQxSEbpl0EG3WSHdXMsqk1/oxwxL7B8et+J1GfXGkkyTev99ulBUljba5d0AWtXpaUkOnv6650ag+SOjaSPpBjs6R/GwcHRhqUYSOtEDdQuiMhKpdR5EVkJ4PWgDZVSzwJ4FgB6e3sDAzpCCCkirZh9BMko/TLIoGAs7L7LbNOfhSFHI0mmaX2QtBSIJoVNuoC1q9JSEiNBckc6OpYHXutEiRqo7QbwUQCfq/5+IbYWpYjJBCSKmQgpBjSFISZanWVqxVHRFGh57fBmgs6OXgaAmrwSqDcw8fZtcxxZ2vRnPfuShVthI0mmaX2QtNS77t57TBJO0/kNmmmMSxLpkhsliYcnHugGgFrepP7AoNGy1yeDtvF/5/0OspTmuoepvIjXP0h82NjzfxMV45CbRWQIwFOoBGjfFpGPA/gJgMeSbGRSmHKQmuUl2eS12RTSZv6Te/CaEA89qMlylinIqt+/jSevBCqFko8OXcLYxGTdTb0rs2VBAVkZZ1/C2vY3+hyThFNHL0wMNC4unLQkkuQX0wMNvY+alrfvHqj1Sb0otr4M1PdJryi2vp1/TCPZwlzDdLBxffxIwJ8ejLktuSBsXluj9xFC3EQPavJQDFo3Ljk6dKnhNo2OI41gLigg4+xLPAS5g9a581l8BtC6JLKMwTdpQhR5JKV1pMS0bCZCCCFFQw9q4ioGnQbru5biwKlFABo7rQXRKJiLy2gkKCBjoeT40fuD7s7XrLhwXE/KWwm+s5bCkviIWhRbf8+L/RdYFJuUEgZqIbHNa7N5HyHETVoNztJ2T/Rmwu7pXt6SE2DQcZ8dvYw//M6PYpHBZR2QFd2Svr7o+PKGLpF6nlES5yLstdaDM8DNWoGkHptcMt0pMkxRbL3vHhwcZVFsEoqi5DkyUAtJlLw2220IIcUg7XwwUwFZIL4b3V19Qzg4OIp7upfnXpZY9PyrIGfIIJdIoL54dpbnQg8yt9zfhS0xFaEnyVH/YMCiH4Uoiq334yRcSUmxCd03HYWBmo9mph9RzUQIaRX2s/yQRl6bf9bOL1WLsyi2XwqaZ5K2pE+LIGlgo6Lj/r7hFS3e3LMCQPaFs/Wbc8od84FXABtAYD+KIn1s1I8JscGmb+YBBmo+mpl+RDUTIaRV2M+yxxTUmNalkdfWbNbOXxQ7SAZpM/uXZZ5es1ylsEFCUZzKWjHq0PsGMC0t27H3BLZtXFM3C2dzTuOapcxaFkvCY5I1+mWMUaSPhLSKTd/MAwzUCCG5RkRuBPAKgHNKqQ+KSDeAnQDaAPQBeFwpFUtkawpqsrK9bzZrF+T+57+pd93VsllAUnQpYxCtGHV479Vz1HbsPVEL1rzzaCsXKsosZRKkOT5ljo07I4tikyzIcV9joOajmelHVDMRQlqF/SyQTwI4DmBx9fXnATytlNopIl8C8HEAX4xjR0EFqP3r0qBZMWt9dqKRDLLRbFnapigmmgUkZQ0SgmafbNwSTe/dtnENduw9gc09K2rn21Yu1MosZQncHVPWLHS7AAAgAElEQVQbn7JAL3AclEvm36ZRUeygQtp5NYKIG11BENYsqmzY9M08IEqp1HbW29urXnnlldT2lxXMY7OH5yX/iEifUqo3o313APgagD8D8AcA/hWAnwJoV0pNich9AD6rlPr1Rp+T9djUakD09L4TeOalV/HJB283Bl3e3wHUDBqA5jdAzT6XzCTrwEO/1mGum/e+e7qX12bRvOV7upfPkETGdRMdtb02ZDk2VfdfiPEpDfSi2ABwR/uimmzNW95yf1fppbHeedLPD89LPrEdnzijlgDMY7OH54W0yF8B+CMAi6qv2wC8qZSaqr4eAnCb6Y0ishXAVgDo7OxMuJmNaVU+GVUG2Uza5ros0kWyLvIcVRJpcg61kUS26qZW8ELnhRifMoESycbwnJQGBmqEkFwiIh8E8LpSqk9E3hv2/UqpZwE8C1SeWMfcvFCYAqIws2zNzD6CZJCP39eJe9/Zhg2r2/D0vhMz9pWnYt+u0CzwSLqOW1RDjiD5YpAk8urk21h9y0I8dteqOklR2BnFohqIFGl8SgPPGdJzIvW7Q7LgdYVNPe04dn4Mm3tWBDpmkmLBQC0BmMdmD88LaYENAB4Skc0A5qGSA/IMgKUiMqv61LoDwLkM22iFKSBKyqREd8LqO/0mFs+fjT39w7VZoGb78geQLuSxuUSzwCNv5id7+odxcHAUYxOTdQWIjw5dwtzZN87YPu4ZRZvA1tEi5oUZn1olbFFswOwOmUfHvrjZf3J63PDOT1KFwItSMDrvMFALCXOq4iWt88frVjyUUp8B8BkAqD6x/kOl1G+LyHcAPIqKs9pHAbyQWSN9hAlq4pAdNtrf2o4lAIBnXnoVazuWYMv9XVb70gPIR9d34A+/86NcBBxRiTvnzCXzk1ABTlVqtbZjSW2WTZdEHjg1gm0b12BsYhKPrrsNi+bNjuXYggJbve16G1zpi3kcn5LiuZcHa8H7gVOV0hBjE5N1DzTGr0zWveeta1Mzlv3blBGTRDmpMcTmupHkYaAWEpucKuaouQfPf6n4NICdIvKnAI4A+HLG7akRZpYsDtmhaX9PPNCNxfOnb6K9HKNf+flbrIIQPYDc1TdUM5soYH4RgPhniFyq42Yzu+f1F38w94UPvWdGoOT9jtMUJCiw9bfdMzzJOvi1wNnxKRUCcqsWzZtd9/qmubNmLPu3KSOZjR/MicsMBmqEkNyjlPoBgB9Ul08BuDvL9gQRdZYsqrzQtL9VyxfUgqxH13fgCx96D557eRBjE5M4O3oZABruS79R0D+/qHKYtMwuspDv2czu6dd75ZL5xjbqQdue/mGcu3gZ23cP1M1ARp2ZDLoxNbV91fIFmQe/JvIyPiWFl38GBOebMUfNPWyuG0keBmohscmpYo6ae/D8ExeIOkvWaCZOD+K8bb2b1qD9+T9v8fzZtddHzlzE0aFLVhKXMpiNpGV20UruWlAQ1Cw4Cvt0PqiN+ufsPzlSZ7PuMX5lErsOT6djtXpOXZqZJI3Zf3KkLr/x4ODojLwqfRtvO/97duw9UVcmooz1w9Kso2Zz3UjyMFALiU1uk2kb5kRlC88/yTONZuL0oAuAlbTS/3ne77GJSRwduhRPo0koWsldC5Jnxi3btGljUCkILyeSlA+9TwQVHva2CSp4bSoTkUX5i6wxPSxJ6jzYXDeSPAzUMoLmFhV4HghpTqOZK1MQ1+yfadDnre9aigOnKiWfgiQudHlMhlZmiILkmXHLNsO2Ue9Pj6y7FVenrgMIZyfuqKMjCYF/VtrUd0wz183KRGxYXb6ZnQ2r23Dg1Egq58HmupHkYaCWETS3qMDzQIpEFkGMP+hqRVp5T/fymtRlT/9w3RNt7yY5qbIBcVOmG/wgeWYWNcr0Wbxj55fXSaeC+pZp2ZNz5a2cAUkGU5mIMtr18zyUDwZqhBASE3kJYvyYLJ8BaDfcY7Wb5DjKBqQBb/CzwS+XMgVhgLlv6cuenMulcgbEAeg+WIHnoTQwUMsImltU4HkgRcLlIKbRbJ9J0nZocKTm+PX4fZ24951ttffqrpGuzlLxBj86rcxGNpJLmfqWLuF6/L5OvGvlYgDT18/GdZIUH71MhB78b989UFs29QnPUMfLfcu7AUnU80DyCwO1jPDnYtnkahUxnyvKeSDEVVxyQfQHZmFn+/b0D9ekNX2n38Ti+dM1jEyf5VruGl0Bg2nmBpn0bKTet54//FpdP2uETbv0IPPGRTffGlujSaaYvs/bdw80LcisF232yHPB5qjngeQXBmqOELWQdtEowzESkgb+YCrqbJ/n1tfss/Iq+ywjzdwg05yNXL1iIe7sXFZ73Wq79GDuxpuWrkys4cQtKAWswPNQOBioEUJIgdBnFACzhMwGT2Kj3wyPTUzi0OBIXa6Rh8uyT1I/i7app72hG2TSs5FPPNBdt6zXffPYsLoN23cPGO3aG6EHc3/7b948H3PTiUNEKaRdxILNLExdbBioOULUQtpFowzHSEiShJnZMskV9XX+93/lh6dr0hqgftbDJdknmYl/Fq1VWVTYPDa/3LKZS6Uu5wJQ1+900xF/n9P74dvjb7wW/QiJ65gKMvuLYo9fmawrpO13SEyzgHRSsDB1sWGg5ghRC2kXjTIcIyFJEmZmyxTUNQ30qtKatR1LOHuWI+KuqRY2j63l4ttav9u2cQ0NYsgMh1FTUexmhdbTLCCdFCxMXWwYqDmM31gDQCkNRwgh9oSZ2WpULPvR9R11s2ubetpx7PwYHr+v02j64JqZCKkn7ppqYfPYggLFIGMTT841+fZ1rF21FI/dtcpK+kgzkfLg79OmotheofUg6WOaBaSTgoWpiw0DNYdpZKxRZsMRQkg8mKz29UDv6X0narNrAHBwcBT3vrMNi+fPrlsPVCRGuw6fo9tYSQhrmx8UKAbNtOlyrt/4pVtxd3fls5956dWa9PHAqZE6mdumnnb8+feO4+jQJewdGKaZSMkwFYPeeWioofSRBaSJ6zBQI4SQEtNI6thoxg2omIt4N9nNJEakuLRi5x8002Zar8/imWRux86P4ejQpcoH0P2uvEjAsu17CHEIBmoOYzLWoOEIIcXABang2dHLGJuYxJb7u5q6/wGYsawXLv6dX+7GqxfeKk1uRLNaZEVHlxh6fShKfkzQTJtpvd4fTTK3zT0rcHXybay+ZSEeu2sV9v3vdH0sMn4n0/Erk1jbsQSPrLsVOw8NAai4IO48NBQofTQVkNZdSUn6FMHgJU4YqDlCXLllzElzD+YNEhMu1B3zZGeffPD2po59elDpvT538XJgUeyi07I5Rs7xz6J98sHbcXd3W0uSsWbBr34D50kcdcnai/0XcHToEubOvjH6gZHcoBez1l1Br05dt5Y+Jl2KgoTHu657B4brrl1ZJfUM1BzBlFsWJUeNuAfzBokJF+qO2bbBH1R6rz25o6kodtGJ20UxbyRRFLtZ8KsHh36Jo+cGqUsimaNWIoLkjpQ05hdeOwAM1AghBNemruPpfSdSlSFGqTsWt1zStg3+gK7RTfq5i5exffdA4aUqcbsoBuGqxDKJmYhmwa/e73QDEd3A5Asfek9t1o0Fr4uNv9CzX+6oLxe54HXR8ByGN/esqF3HMl83BmqOYMoti5KjRtyDeYPu88bPruKZl1513rHQk4Sk3U5/QNfoJn3X4XO1ZZfPZV7Q5V1Asc9ps+A3qN8FSdlY8LrY+As9N1sGKtJHPcg3LfsfiDBnKl32n5w2JfKuXZkLeDNQcwRT3pJNLhNrrbkPz29+OHLmIs6OXgaAzI0+TIxfmaz77TLjVyaxffcAgPqbGldniFxFv9Z5uO6EpIW/0HOjAGz8yiQWzZtdKyMC1Oe16ctA/QMR5kyli2nmvIzycg8GajmHtdZImRGRpQD+E4BfBKAAPAHgxwC+BaALwGkAjymlLjb6nJsXzsUd3ctxcHAUu/oqUgsX860WzZtd99s1/G5ppnyjsptwhEW/1q5ed2ImrvGJmGlU6Dlo2Xt4BCB8LhtzplKBBi/13NBsAxFZJSLfF5FjIjIgIp+srl8uIvtE5NXq72XJN5cQQup4BsD/rZT6BQDvAXAcwJMAXlJK3Q7gperrhsyZdQO2bVyDe7qXY8PqNmxY3VZbdolNPe24p3u5s3p978bpqYfejU097bijfRHuaF9Udx43rG4zridm9PPo6nUngcQyPpH40L9PH76rw7js/555426jbQhJCpsZtSkA25RSh0VkEYA+EdkH4GOoDDafE5EnURlsPp1cU4kJ1lojZUVElgD4ZVTGIiilrgG4JiIPA3hvdbOvAfgBLMYmTxe/p38YR85cxNGhS9jTP4yVS+Y7I4P02pgHvb6eP6LnhYxfmaxb7x0HJZFm9POYh+tOKsQ9PpHm2OSS7ekfNlr368t+C3/mTJGw+Gv86YZHYWkaqCmlzgM4X10eF5HjAG4DwMHGAaLkPzFnihSEbgA/BfAVEXkPgD4AnwSwojpuAcAwgBWmN4vIVgBbAaCzs7OmgR+bmJy2/oYb9c48XLD0t0XPHwGmZZCelb8fSiLNlL0MQI6JdXwizfHX9gOajCWW0kfmTJGw6P/Pjp0fq/XHKP/XQuWoiUgXgDsBHITlYEOiY2P6EcVMhJAwOGw+MwvAOgC/p5Q6KCLPwCcjUkopEVGmNyulngXwLAD09vYqTxd/aHCkZvnsyVsOnBpxQqYXxdI/K/T8Ef2cPrLuVlydug6g3nJ5w+o27B0Yri2TCmmVASCxE+v4lHRji8CG1W04cGoEm3tW1B5u+McSG0t/v6yROVMkLPr/s809K2rromAdqInIQgB/C+D3lVJjItOPHxoNNnwqFB0b0w8WxSZJ47D5zBCAIaXUwerrXajcCF0QkZVKqfMishLA6zYf5slmzl28XCeBWTx/dm7khnHXWWu1LZ70Q5c7BkmMKPGLH8pJMyXW8Yk0Z0//MA4OjmJswiyv9l6HlT4SEhZTP4var6wCNRGZjUqQ9g2l1Herq60GGz4VIoQkgVJqWETOisjPK6V+DOBBAMeqPx8F8Lnq7xdsPs+TzfileXmSG7ok09SlH3XnNEBiRIlf/FBOmh1xj08kBLbujGFdHwmJQot9q2mgJpWpsy8DOK6U+kvtT7vBwSZRbEw/WBSbJI3j5jO/B+AbIjIHwCkAW1Bxs/22iHwcwE8APGbzQbpsRpfm5Ulu6B2DC9JBXfqhyx2DJEaU+MUP5aSZE9v4RJrzxAPdWDx/9gwzEf82Ho1qr3k2/o1qsjUqkk3Kjb+feWYiURClGk9yicgDAP5fAP0ArldX/zEqeWrfBtCJ6mCjlBpt9Fm9vb3qlVdeidRQEg1TfhHz2kiciEifUqo363a0Qm9vr/rtv/gbPPPSq7inWk8NALbc31Xbxht4XZEWmnh63wk889Kr+OSDt2ceXHptAVB3ToPOb9BND+V70dGvgQt9Im2KMDYBvHdKm+27B2oz0Xe0L6pJ2PRl/+st93fxQRMJhe34ZOP6+DKCJ+4eDNswki6m/CLmtREyE5OzF1BftBlwsxC2h0syTV3KGPTkGoDRGYsFsuOBclJCWqSRPJJySZICoVwfCSGkqJicvXS3Qk865oq00IRLMk2/lFFPojad3yCnNsr3okM5KSHhaeYMefnaFBbMmdXQKZKQuGCgVnBM+UXMayMkGN05UXductkB0iW3RxNB8sUgBza6QUaHUlFCgrEpRGzjDAkAL/Zf4NjUBJsi5KQxDNQKjinfjDlohASjOyf6cUlaqOOS26MJK/ki3SBjgVJRQoIJXYg4QPq4tmMJtm1cwwLYTQhdhJzMgIFazmlWjNjGTISBmz08d8XH75yoF792SVqo45Lbo4kg+WKQxIhukNGhVJSQYGwKEdsUxf7j37gDd3e3cRatCTZFyEljGKjlnGbFiMOaiZDG8NwVn/0nR2ryRgB1cjzdEMMkmckKvc0u3jgEyRf19bqMyH+u/ZbYLp1716BUlJBgbAoRsyh2fNgUISeNYaBGCCEafnmj9xQQMDsUujDD5qok0yNIvhjkDAmYz7X/tQvn3jUoFSXEEhvXRhbFjgeeu8gwUMs5zYoR25qJEDt47oqPX97oye78DoWAOxIOVyWZHkHyxSBnSJMbpOe05tq5dw1KRQkJxqYQsV/6+Pzh17D6loV47K5Vtdl9oDJOBRlk0ESjgk0RctIYBmo5J0qOlE0emx/mZlUo63GXlSAHyOcPv4ajQ5dyIy1z1RUyrBskUJFIuizzzCv6jSWlpaSoNCob4mEaf+7sXFbLSfMKYh84NVI3Numf+9zLg/jKD09j78Bw4DZlwFT2hoSDgVrBscmpimsbQoqG9892bGKy9hQVAFYunYdf+flbciMt04/DpRsFr10eXtv0c/3Wtana8h3ti7Bt4xrjU3DSGv4bS9f6CiFpYRp/9HXesj426X+33YYQGxioEUKIBYvmza4tty+e77TUMO/o5/qmudP/pu59Zxud1pKGuSSk5JjGH32dt6yPTfrfbbchxAYGagXHJqcqrm0IKRqevv7R9R04f2kCR85cxOpbFmJTTzue3neiNqvjoqxQRz8Ol9BzFfy5I16OyCd+9efQd/rNurwQXZpHyV486LkknLEkZaZZHpsp72pTTzu27x6oLY9fmcTajiV4ZN2txrIjAPPYiB0M1AqOTU6VTU6aDcxjI0VD19c/9/Igjg5dwp2dy7Cnf7gmJQTgpKxQx1WzkSDjCz1HpO/0m3jqoXfPyAvxzjcle/HAXBJCKjTLYzN9V7zxCUBd7trVqeuBtvTMYyM2MFAjRhrlpDGPjRCSGUHSPEr2CCEuENbSn2MXaQADNUIIsSBIPri+aymOnV+OTT3tzror5g2TJNIkNzo0OFInMXqx/0JNXkRJJCEkLfyW/p7cUV/2Sx839bTj2PkxbO5ZEbgNIQzUiBFTThrz2EiZMckHK1K8RTg+PF4LHlyXQeYBkyRSP/+e9NQkMfLkRZREEkKSwJRbNn5l0lhSRF/2Sx/3nxzBwcFRAKhtU4ayI3pZFj0HkA/SzDBQI0biqM9GSCmgbCU7bCRGvD6EkBjZ1TeEZ156FQdOTQdaazuWTG9gKX301Bl6wFcGE59dfUO1fL5j58dq59DFPGoXYKBGjCYgJjMRGoXMhAYq5eTs6GWj5G7lkvlOuisWQZKpP8U2Oap9+K6OOukjXQzDQakoIXZsWN2GA6dGsLlnRc1Q6pF1t+Lq1HUA9tLHshr4bFjdhr0DlcB0c8+K2jpihoEaMZqARDETKSM0UCknz708iF2HzwGol9w99dC7nXwq6GrB6zCYHNJ0RzVPYuTJi8p6ExQVSkUJsWNP/zAODo5ibKK53LGR9LGs6K6+/nGbzOSGrBtACCG5htK6dAmSFfE6xAPPIyF22IxF/D41huenKZxRI0YTkChmImWEBirZIiKfAvA7ABSAfgBbAKwEsBNAG4A+AI8rpWKd7vQXRNVdCPefHHFOPuZqweswBBWZNS0fGhypW88C2c2hVDR+shqfSLKEGYv8y3pRbNM2cRS8dr2QdrOC4qQeUUqltrPe3l71yiuvpLY/khw2eW1x5WsxD8xtRKRPKdWbwX5vA/AygHcppSZE5NsAXgSwGcB3lVI7ReRLAH6klPpio89qdWzyip3e0V5xgPR+b7m/i/KxFPFfBwAzroV/G16j4pLV2FTdtzPjE3EDvSi2aYwCEMt4ZBoHOc65h+34ROkjiUTYvLa490VIlVkA5ovILAALAJwH8D4Au6p//xqAf51aa8T3m2QD3SCJG7g1PhF3SEMeyTGuEFD6SAjJJUqpcyLyBQBnAEwA2IuKlOhNpdRUdbMhALclsX/dSVEvXPpi/wU8fl8nXr3wVq7lHC47RQbJF00FZOkG6R56HSWXJFlxkvX4RNwjSlHsqPthIe3iwECNRMI2ry2pfREiIssAPAygG8CbAL4D4AMh3r8VwFYA6OzsDL1/3UkRQJ0L2LtWLs69zMRlp8ggh8JGbmx0g3QH7/p5uNa/4iDr8Ymki00RZ5PboX85DvfDIhbSLnORbAZqJBKmPLGkcseYk0YC+DUAg0qpnwKAiHwXwAYAS0VkVvWpdQeAc6Y3K6WeBfAsUMkBiaVFlJqki42skdeEZIN74xNJjNBFnBMco4pYSLvMRbIZqBEjzQw80jQTKRo8T7FxBsC9IrIAFWnRgwBeAfB9AI+i4qz2UQAvJLFzv5Oi7gI2fmUS23cPJCbrSkOW6LJTZJB8sZkbWyM3yDxQFMmg7vqmLxeMTMcnki42RZzTkj4WUTVQ5iLZDNSIkWYGHmmaiRQNnqd4UEodFJFdAA4DmAJwBJUn0N8DsFNE/rS67stJ7F//ZwhMP9nb0z9cK4YNJCPrSkOW6D8+lwi6EQlav6d/GF/54WkcODVSkwJ5yy5KO4MoimRw1fIFuW27LVmPTyRdbIo4pyV9LCJlLpLNQI0QkluUUk8BeMq3+hSAuzNoDnEdSiJJinB8Kik2YwvHouiU7HwxUCNGmhl4pGkmUjR4noqH3wHSk7ck5bRlI0t02bUxKVpxg8wDcUkGiyKhJMQVbIo4+7fxy7PHr1SMqfwS7bBybdcLXkehzEWyGagRI83ypmzMROLKYytaTlfe209m4neATFrGYiNLdNm1MSlacYPMA3FJBosioSTEFfzfTdOY0mibu7vbaoWqTRJtfbnZmO4fBz3y/D23Ob9FhQWvSWLElcfGnC6SV8avTOLpfSdwdvRyrJ97dvRyIp9bGEK4QR44NYLtuwdycS7Pjl7G9t0DuWkvIcQeb0btrWtTtXWmZW+7MJ/T7D3EXTijRgghLdJIivjMS68CiNdGeFffUNPPddm1MSnCukHu2HsCBwdHcXx4HIvnz3bWPMVDt6hupb0lcV0kJFcsmjcbAHDT3Olbc9Oyt12Yz2n2HuIuDNRIYsSVx8acLuI6fimiJ9E4NDiCY+fHYrcR3rC6DQdOjTT8XJddG5MirBvkto1r8NQLAwDyYfWsW1S30t4yuC4WiQtjV1rOWyLuE5RLGzav1vRgig9k8kvqgZo/3whAofKPSPywT5C8sv/kCA4OjmLH3hP4wofeAwCxGHz4P5c3aNHYf3KkLp/QlHjvkvGGTXtJ8Xh9/CqeemGgpbwl4j6Ncmn15WZ5tUWso1ZmUg/UGuUbMf+oWLDWGik7j67vwIFTlaBqV1/lqWgcUkj/55Zt5iwuHl3fUTOAAWCUFcYlN4wDm/aSghJk514yq/JSwGtNNJoGaiIyD8A/AJhb3X6XUuopEekGsBNAG4A+AI8rpXjnTQghVVYtX4BtG9dgx94TNalaM8li1M8l4dElgIcGR2plFfRzGpfcMA5s2kuKx01zZs2QwOW5zAQxE5f0kRQLmxm1qwDep5T6mYjMBvCyiOwB8AcAnlZK7RSRLwH4OIAvNvswU74R84+KCWutETItU9zTP4wjZy7i6NAl7Okfxsol81uSQXqfu//kCOUtMRAkKxy/Mmkst5C1JFJvL/tAsXnr2hRe7L9Qu97+5SKMA973afzKJBbNmz3jO1XE2mB+vDEdQMNrvWPvCWzbuIY5iyWhaaCmlFIAflZ9Obv6owC8D8BvVdd/DcBnYRGoMd+oPNjUWiOk6HjOg2MTkzg6dKm23sa50eZzy+TqmCRBssK1HUuM22ctidTbyz5QbG5ZNNfqxjzP6N8nYOZ3yhsvPcm3aZu8411DvzutyanW+w0Ax86PGZeLdG7KjFWOmojciIq8cTWA/wDgXwC8qZTyijQMAbgtkRbmEJsCzUUr4hwVmsuQouMldutSNU+60ooMsoyujkkSJCt8ZN2tuDp1HQDqJEdZSyLp3FgeViyeh7u722YUSDYt5xXv+3T52hQWzJk14zvlOd1u7llRe0BRNMlvkAmI36l2x94TdefBvwwU79yUGatATSn1NoC1IrIUwPMAfsF2ByKyFcBWAOjs7IzSxtxhY5hBU40KNJchRceT7Jy7eLlOQrd4/mwnJEte+8ouldGljLrcUXdd06WPlB5Oo8vSKLsiUdjTP1z7Pnmv9e+UyRGxmfthEYnLGZLkh1Cuj0qpN0Xk+wDuA7BURGZVZ9U6AJwLeM+zAJ4FgN7eXtViewkhJFd4kh2/hM4V6WKrEsyioEuv6q5VgNMapYfTmGRpZe5LJEHofFiBzpClwcb18R0AJqtB2nwA7wfweQDfB/AoKs6PHwXwQpINzRM2hhk01ahAcxlSdHTJji6hc0W6aFM8uwzoUkZd7qi7runSR0oPp9H7uPeakDB4BZl1MxH/3/1FnDf1tGP77oHachlmc7N0htRVB2U53y5gM6O2EsDXqnlqNwD4tlLq70TkGICdIvKnAI4A+HKC7cwVNnlV/m2Ys5Y+POckDUxOXv5ixUA8hbBbaV/Z5Xu6lDHIYU+/bnqSv6sFstNC7+PsSyQKzR58mPK3tu8eqM2Cl6Xwd5bSx+deHizd+XYBG9fHfwRwp2H9KQB3J9GoMlLWnLUsc9TKes5JupicvID6YsVAPIWwW2lf2eV7upTR77Rmum66u5qrBbLTQu/jRXAgJDmkbJK/rKWPZTvfGRIqR40QQkg4TE+CTcWKs5IfuiLBzBr/E32T65p+3YLc57J2g8yCILc6QpJkU0977ftYlmLQWUofy3i+XYCBmiOUNWctyxy1sp5zkg26u6LucOaSA2QUiugaGSRf1K8b3SDdga6T5cQvV3Z9DG2W42VT1DtqUew4ZNgune8yFED3YKDmCGXNj8ryuMt6zkk26O6KfvIsPyyia6SVfJFukM5A18ly4pcruy671ccVU2Fqm6LeUYtixyHDdul8l6EAugcDNRKJNIt6l8H0owzHWHb87op68es8yw+L6BoZJF/0S3/oBukGdJ0sJ43kyi6ijyumvmpT1DtqUew4vhMune8yFED3YKBGIpFmUe8ymH6U4RijIjATlz4AACAASURBVCLPAfgggNeVUr9YXbccwLcAdAE4DeAxpdRFEREAzwDYDOAygI8ppQ5n0W4/ursiYHaAzKPlcRFdI4Pki0HOkCY3SM9mPI/XNG9k6TpZlPGJJI9JOq3LpuMq6l2G4uBlOEYPBmqEENf5KoC/BvB1bd2TAF5SSn1ORJ6svv40gE0Abq/+3APgi9XfmeOXN3pPAQGzk2BeZtjyLNsMIki+GOQMCZivof91Xq5p3sjYdfKrKMD4RFKmkWtiXI6KZXBmLMExMlAjkUizqHcZTD/KcIxRUUr9g4h0+VY/DOC91eWvAfgBKjdCDwP4ulJKATggIktFZKVS6nw6rQ3GL2/0JCR+J0EgXxKOPMs2gwiSLwZJf0xukJevTWHBnFm5vKZ5I0vXybyOT2UyY3AFvYi3KcfLVNTbX/jbdj/654xfqTxcOjQ4Enit4+oPaRXFjutc5YHCBGrM8UmXKEW949pXWvlxafYp9tfQrNBuboYBrKgu3wbgrLbdUHVd5oGaTpAD5POHX8PRoUuJSreK6NKYFmHdIIHs3dFahY6KkXB+fPKKF+8dGK7rr8yvTI5mOV5xPXDwf45XGNwrUu2htyWu/pBWUewylQQpTKDGHJ/ykFZ+HPtUPlBKKRFRYd8nIlsBbAWAzs7O2NvVCO+f2djEZO1pJwCsXDoPv/LztyQq3dL3zZuycOg3IcD0TYx+Dd+6NlVbvqN9EbZtXOO8G10j/Ddw7DfhcHV88vqs3l/1fkyKg821jqs/mMZC9qvWuCHrBhBCSAQuiMhKAKj+fr26/hyAVdp2HdV1M1BKPauU6lVK9b7jHe9ItLGNWDRvdm25ffF8fOr9azhjkTP0a3jT3Onnn/e+sw13d7cV45qWIBckRpwfn7w+q/dXvR+T4mBzrePqD6axkP2qNQozo8Ycn/KQVn4c+5TT7AbwUQCfq/5+QVv/uyKyE5Uk/Usu5Kf58fT1j67vwPlLEzhy5iJW37IQm3ra8fS+E7XZlyQkivq+STj0HAh/vomXo/aJX/059J1+sy4vRJcM5k1KqOeC5HlmMGUyH588ma7nPurPNypTjk/ZMV3rTT3t2L57oLY8fmUSazuW4JF1txrLjoTZl0eStdbKlGOZi0DNJleIOT7lIUrOWhz7SRLmWAYjIt9EJTH/ZhEZAvAUKjdA3xaRjwP4CYDHqpu/iIr19UlU7K+3pN5gC3R9/XMvD+Lo0CXc2bkMe/qHa7JEAIlIFIto/pEWQSYjeo5a3+k38dRD756RF+Jdx7xJCcuUCxIFV8cnv0wXqM834nUtD6Zr7Y1PAOpy165OXW/J7j6tWmtlyrHMRaDGXCHSCFP/yFufyVt700Qp9ZGAPz1o2FYB+ESyLSIkBEGSQUoJCwHHJ5J7JGA5D+StvRHIRaBGCCFFJUiKuL5rKY6dX45NPe10anQckyTSJDc6NDhSJzF6sf9CTV6UN0kkcRtPjuuViYgiYyPFRZdrf/iujprcUV92uc9s6mnHsfNj2NyzIhftbYVcBGrMFSKNMPWPvPWZvLWXxIdJiliRzC3C8eHx2k0+nRrdxSSJ1K+rJ2k1SYw8eVHeJJHEbfafrLdiz2t5CNI6pnyu8SuTxpIi+nIU6WNa7D85goODowBQa28cfdzF3LdcBGrM1yGNKELOYt7aS1KgBJKO0mEjMeJ1JzHw6PqOWvmPRfNoHlRmdvUN4ZmXXsWBU9PBzdqOJdMb5FD66PVnPaCKo4+bztXi+bMzzevORaBGiI7feAPADCMOmnNU4HnIF2dHLxulcSuXzI/VqZFSyuTRn8yaHNU+fFdHnfQxD+6KlGfmhyDjG1I+Nqxuw4FTI9jcs6JmVPXIultxdeo6gHxKH5MywzGdqw2rs51VZKBGckcj4428mokkBc9Dvnju5UHsOlwpq6RL45566N2xPtFj0evkMbmS6Y5qnsTIkxflwYWP8kxC8see/mEcHBzF2ERzuWNepI9JYTpXWZ8HFrwmhBAXyYkEhTQhSFaU5+ub57YTUlZsxiJ+tys4dB44o0Zyh8l4I+9mIknB85Av/MVCdbfAIKlZFBkji14nT1CRWdPyocGRuvVhZYVecWNvv0nJEfMgzyTFQO/TUb4TZJowY5F/WS+Kvad/uJbzWJRr4u9nALDl/q6685B1IXiplPVIh97eXvXKK6+ktj+XYK5QtpjOf9GuSVbHIyJ9SqneVHaWEK6OTV5R0jvaKw6QW+7vmiE187Yx/Y24jf/6Amh4rZt9DgD2A40ijE2Au+NTkuh9Osp3grSO6Rp4FOWaZNnPbMcnSh9TgrlC2VKEotjNKNrxEA2HZBgkASg9IiQYfieyx38NinhNHD0mSh8JIcQhdCmjXtRTdwjUiUPGSBfIZAlySzQVbfW7QdpgKridFWnJMEmx8RdkDvudIK1jKortFVAvyjXJQz9joJYSzBXKliIUxW5G0Y6nrOiOjADqHKhM7lOmgtmt7DPPMhZXCXJLbOTGFsZpzCUrdu9YPVxpF8kXesHuF/sv4ODgKIt2B5BUPt+e/mGjG6T+OglHxDSLTuehnzFQS4m85z/lnSIUxW5G0Y6HaDgqySAhsSlyzWtNSK1gN1AffJCZ7Oobqj0cOXZ+rFaoOdYizSlKH9MsOp2HfsZAjThFUoYYWZqJFM20hCSLX8qou3WNX5nE9t0DdU8X45At0gUyWYLcEpu5sXlukJ7TWtBTZZfkhi7JMEl+8c8SuzTD4RobVrdh70BlzNjcs6K2rlWykj6mWXQ6D/2MgRpxiqQMMbI0E6HJBwmDX8roLe/pH64VwwamJWVxyBbjkE+SYIKKWQet39M/jK/88DQOnBqpkxsBZimhS3JDl2SYhJQBk0QxDkliVtJHF4tOZwkDNUJI6bk2db1WLybrGYkgxq9MGpddwp8rkXR+QVHxru9b16aM6wlJgyi5QmnmFxEDSUkSA6SPR85cnFEHspWakA33WVIYqBGnSMoQI0szEZp8uM/Fy9dqMxJJaeFbZdG82cZll2SLQbkSrp5TV/Gu701zZxnX+6HckCRBlFyhNPOLSAX9Ox9nnpX/c/0Fr3fsPYGDg6O130D9uB81X84kCS/zuFaqQC1KrhDzi9IlqfNrYybiv9YAYrn27DPus3DuLHS1LwKQnBa+VfR8AT0vwCXZoj9XIun8gqLit+338kKC8kEoNyRJECVXKM38IlIhqTyrZp+7beMa7Nh7ou5a+5eB8Nc/SBJeVkoVqEXJFWJ+UXlodK157YvNz65O1bTwrVrz6uYeAGKrT6bbCO/pH3ZSVuS3Oo7rnLZKUB0zl0w4dPafnJ6N0PNC9OtOaSlJmii5QswvKg+Nyovoy7z+rVGqQI0QQkwsWzAHv3l/FwC0LBnxpD8e3nKrs166jTAAJ6WafqtjL5DIWpZpkmN96v1r6qSarp1HAHVun57skdJSkjpRcoWYX1QegsqLsA/EQqkCtSi5QswvKg+ma81rXw7mzLohNumYJ/3x5B76civoMpRDgyM1GaRLsiJXrY51OZb32vvtSTVdO48m6Y9+3SktI0kTJVeI+UXukZTBi1+iDUzb+HvLSVn4RyGvRjelCtSi5Aoxv8g9mDdYgefBTTzZ2v6TlRtqL9n6Cx96T2z/EHTb5EaykjhqrBUBXUroXZu7u9tyISfVCbLLjiotCpKEEgJEyxVifpF7eOU79g4M10mpW304mTfpY1LnIWlKFaiRYpBmrbUk9hMXZcmfFJHnAHwQwOtKqV+srvv3AP4VgGsA/gXAFqXUm9W/fQbAxwG8DeB/VUr9fZrt9WRr3m9Pbrerbyh1eZouwyyzNE6XEuqOaHmQkwYSg6woSBJK7Mnb+ERKTBrW/XmQPrrargCsAzURuRHAKwDOKaU+KCLdAHYCaAPQB+BxpVRx7xYJIVnxVQB/DeDr2rp9AD6jlJoSkc8D+AyAT4vIuwB8GMC7AdwK4L+KyBql1NtpNdbvwug5Y8UpTwtygPTjl2GWlUYFp12Xk+ro112XGEWVFgVJQkkovoocjU+kfJgkinHIEfMmfUzqPCRNmBm1TwI4DmBx9fXnATytlNopIl9C5QnRF2NuHyEzSLPWmss5amXJn1RK/YOIdPnW7dVeHgDwaHX5YQA7lVJXAQyKyEkAdwP4/1JoqhFdChlW/hHkIOmX7AXJ1sLuu8xSSf2cpuVS2ciN0rQ+yFUzqmwzSBIa5/G56KoZJ3kfn0hz8prb5GFykY3ju276XH1cerH/Qk36v23jmsCi2P76bEnJsJM6D0ljFaiJSAeA3wDwZwD+QEQEwPsA/FZ1k68B+CwYqJEUSLPWmsvkrb0J8gSAb1WXb0PlxshjqLpuBiKyFcBWAOjs7EyscX4pZBiCHCT9kr0geWPYfZdZKqmf07RcKhu5UZrWB7lqAtFkm0GS0Lhw1VUzZZwen0hz8l7E2+8iq6+L+3P9QViYotj+13Gf36TOQ9LYzqj9FYA/ArCo+roNwJtKqanq68DBhrgHTSiiY1MUm+c3PUTkTwBMAfhG2PcqpZ4F8CwA9Pb2qpibVqOVgtRBDpJ+yd6x82NG2VrYfZdZKplF0ehGbpSm9UGumlFlm0kbP7jqqpkWeRifSHPyXsQ7qe950OfqyzZFsccmJnH52hQWzJmVqAw7r0Y3TQM1EfGSZPtE5L1hd8CnQu5RFhOKJLAxHOH5TQcR+RgqSfwPKqW8G5lzAFZpm3VU10UmSzmgyUHSL9XwtvGcJQFzkW2b42hFptkqZZDJ+WnkRhlGkmjjXhnm/MblBpmFnNQV0hqfSPKwiHd0bJ0hPVxzinQBmxm1DQAeEpHNAOahkqP2DIClIjKrOqsWONjwqRAhJG5E5AOozPL/ilLqsvan3QD+RkT+EpVk/dsBHGplX1nKAU3SRb9U49H1HXXOkoBZCmlzHK3INFuljDK5Rm6UpvWNPqeZe2WY8xuXG2QWclIXSHN8IimSM7dAp7B1huQ5nkHTQE0p9RlUHItQnVH7Q6XUb4vId1BJkN0J4KMAXkiwnSRGymJCkQQ2hiM8v/EiIt8E8F4AN4vIEICnUBmT5gLYV0mZxQGl1P+slBoQkW8DOIaK5OgTrTqqZSkH9EsXTTfLq5YvmOEsaWqvzXG0ItNslTLK5Bq5UYaR6Ni4V4Y5v3G5QWYhJ02brMcnkjxRinjn3YAkLmycIZ8//BpWLpmH9iXz68xEdHRFQJKGIy7SSh21TwPYKSJ/CuAIgC/H0ySSNFFypph3ZU+zc8NzGQ6l1EcMqwPHG6XUn6FifBQZXSaYpRzQFq+Ne/qHceTMRRwdujRDOmJzHGnIPIMkeGWWycVJkAxy/Iq9bCsuN8gyyFmzGJ9IukTJbcprceW4sZU+3tm5rHZuTOfYO59A5UHk8eFxjE1MluJ8hgrUlFI/APCD6vIpVGxlSQlg3lWFOIpi81y6jy4TzFIOaIvXtrGJSRwdutRwm0bHkYbMM0iCV1aZXNwEySDXdiwJ9RlA626QZZSzElIHpXwV4iyKXbJz2sqMGiGEFBJdJpilHNAWr4267M1fyNPmONKQeQZJ8Mogk0uDIBnkI+tuxdWp6wCaF3mNyx2tjHJWQoD8FleOm7iKYm/qaa+NZa4V0k4aBmrECuZdVYijKDbPpfvkQe6o40kWz1283JIrWaPjjksWSYljsuhyQ13uqMuN0nJUa+Val0E2SYpLXosrx41J+vj84ddmFMX+8+8dx52dywLzz/SxxHtPWuczSn5cnDl1DNSIFcyjqhDHeeC5dJ88yB11PMliGHmbiUbHHZcskhLHZNHlhnX9IQO5UCvXmrJJkmfyWlw5MbTxZ/WKhbizcxmA+qLYR4cuBRa81seSIMORpNDHItuC3FHeEwQDNUJCYGME4rJZSJS2uXw8SZEHuaOO7tJnK28z0ei445JFUuKYLLrcUJc76nKjtCRDrVxryiZJnslrceW4MUkfH+tdNaMo9lMvDABAoNOsfyxJ85zqY5GtE26U9wTBQI2QENgYgbhsFhKlbS4fj6ukXSTbJLMJKnwchbOjl2tPPfMu34mrmLOr+CVCpmWvb4xfmcSiebOdsbvW5UJAuSVjhBQBm2Lhe/qHZ0i0XSp4HaV9cR4TAzVCCIkZz0o4Lftgk8xm/Mokdh0+V9umlXY89/IgDg6OYm3HktzLd/y22UWzePZLhLz+oC/7+4Yrdte6Bfej627Dlvu7Kss573MkO7zg33sowXzH6LSUd2Urva5ud+TMRRwaHHGvBp3hOEw188avTDZ8TxgYqBESAhsjEJfNQqK0zeXjIRVMMpvtuwdi38+dncvc+GcZBwW1eG4kEQrsGw6ei0XzZhcqgCbZoOcKAcx3bIUoeVc2xcL113rOmvcbyPa6+dvnz4/z8rcPnJpWtmy5v6v2oKnVnDoGajmijLlCSRLlfEbZxqXrxpy0dPD+OSU5E9BMXqnbGTfKSbKRaaZxPGmh3zikmZAeN624Inp94/K1KSyYM6tlu+u45KT6DZH/Zo6QKHi5Ql5fZ75jdKLkXdnk6vkfLm3buAY79p7A5p4VNXVAltetWX6cniOuqxlMD8miwEAtRzBXKF7SOp95u255a6+LpGFG0kxeqWvkG2njbWSaeTNXaURRkvx1mSAQTtqq9w2g9RyKuOSkNJohcePv6y7lPuWNtHLJbPLaXCLp9t4Qy6cQQgghJJ/EJX10UEJJCEmANL7reRtPEmovZ9RyBHOF4iWt85m365a39qaJSSaYtsOjRzM5ol9X7+Ul+WVyrssam0n8iu7iGEQrMkFve9310csf0RP4bQ0DiiInJcXD39ddktSaTCicMc0w0CxXy5ZmpiSmvDb9f1jQuBRnkekw2OThtYIopWL7sGb09vaqV155JbX9kXQpY25TUnXJ8nQuRaRPKdWbdTtawXZs2r57AF/54Wlsub+rJtEyrXMNr40AnG6niWZt9/5+R/siHB8ez93xuYL/PALI/TktwtgE8N6pqJi+c3n9roVBH9Ntxxib90T53CyxHZ8ofSSxUcbcpqTqkpXxXBLSEnmTybiKBCwTQpKhzN+zKMdu854CnVNKHwkhxBKTTNAF6WBcDpBhPzcNmkn8KLurJ6ob5Kaedhw7P4bNPSuw89AQALTsBplkewnJO6bvXJLfNVfQ/x/ZjjE274nyuXmAgRqJjTLmNiVVl6yM5zIPmNwPXXBEjMsBMuznpkEzJ8CiuDjGRVQ3SJNzWZLObh6tuFcSkmf2n5yuu+V95/afHGn5u2aT+5ZVPhdQOUbveF/sv4CDg6NNj9v0nh17T2DbxjV1RabTHLvSgoEaiQ2X86iSIsoxx1WvjRBCYofSR0JSwZv91wOqOBQBpgLM/oLRUYpXx8Wj6zvq6o3ZKCH87zEVxV7bsWT6DQUauxioEedpZqyRJ+ONIPzHACD3x0TsaVViGNUBstmTVBdknXkjaylfVDfIIKe1IDfIsM6QcbeXkLyTlBrAVIDZXzA6SvHquGhWQNrmPaai2I+suxVXp64DoPSRkFRpZqxRBOONRseQ12Mi9rQqMWwmv9T/yenOWAdOjTQsUOyCrDNvZC3li1o0OuimcU//ML7yw9O1vgLAuMwi14S4gU0B5rSKVydFI6m2vpynYwqCro+EEFJWCiQPIQkTJIlkHyLETWy/m3n+DpdgXOKMGnGeZsYaRTDeMB1D3o8pLkTkOQAfBPC6UuoXfX/bBuALAN6hlHpDRATAMwA2A7gM4GNKqcNptzksJolhUo6LujPWJ3715/DqhbewYXUbnt53ojTFopOkmZQvbwW6g9wgnz/8GlbfshCP3bWqLscka+ln2pRhfCL5wqYAc1zFq9PCb34yfmUSazuW4JF1t9aNS0k41mZpvAIwUCM5oFl+Vt7yt7LMqctpPt9XAfw1gK/rK0VkFYCNAM5oqzcBuL36cw+AL1Z/O41JYpiU46Iueek7/SaeeujdNTlklu6ORaGZlM+7rnsHhluSDKZFI4nRnZ3LcHd3W520KGvpZwZ8FQUfn0i+sMl9i5InliX6uKJLr69OXU9c+mjad5rjNqWPhKSMKR/Nvy6pvLs85vMppf4BwKjhT08D+CMASlv3MICvqwoHACwVkZUpNJOQcORNmlNQWVGrcHwiJGVs5I5JjVEZjH2cUSOE5A4ReRjAOaXUjypqohq3ATirvR6qrjufYvNiISnHxU097Thy5iJW37IQm3ra8fS+E9jU0269L78k04Wi2HkibwW6demjLnfU3SB1KVDcLo42UlHX5KRlGJ8ISRN/MesguWMS0sesC2kzUCMkZUz5aGnlqBUhn09EFgD4Y1RkRa18zlYAWwGgs7MzhpbFS1KOi3v6h3F06BLu7FxWc/QLI+PwSzJdKIqdJ/JWoNsvffTkjn43SO/6x+3iaCMVdUlOWpbxqRVsCjITouMveG2zfHBwFH/+veO4s3PZjJIiYUqNRCnQHScM1AhJGVNeWFq5YjnJSWvGzwHoBuA9re4AcFhE7gZwDsAqbduO6roZKKWeBfAsAPT29irTNoSQKkGSn7SkQDb7cUOSyfGpCTYFmQnR8Re8tgm6vGLYR4cu1RX1tlnW+2KUAt1xwkAtBnJq0EByRrOi2P7XRe2HSql+ALd4r0XkNIDeqqvabgC/KyI7UUnSv6SUckZWlLZM0LS/IEnl2dHLVm3S33929DIAYMv9XYUtVlxkF0MbyWCQVNPkLOeXQcZB0P71tgOVPuiCe12ex6e0sCnITIhOI/OToOVtG9fgqRcGAKCurzVaBmb2xayNVxioxUAeDRpI/ghTFLtI/VBEvgngvQBuFpEhAE8ppb4csPmLqFhfn0TF/npLKo20JG2ZoGl/Jkml52hl0yb9/Z5b5Jb7uwoVwOgU2cXQRjIYJNXU1wfJIOMgaP/+tm+5v2uGA2UaFGl8SgubgsyEtIqpqLfNsmt9kYEaIcRplFIfafL3Lm1ZAfhE0m1Km7hn4vTP8zhy5mJthozmIOmTqSGGpWSwaRurn3Pg1Ai27x6om4FMbGYyY7kjx6cWcEOqSmIk65pjOuNXJmvLb12bMi5fe/v69BtCjINpHiMDtRgogkEDcR8bwxH2Q7eJ6uTo5XQACJXHEbQ//fOeeKC7ps3f1VdxzLLdV1LOlC4Rt4thEKa8naRzdsI6UAa1Uf8cLy/k+PB4Xd7Rrr6h2sxkHPlIeXPPJNPYFGQm+UT/ngflfKXFonmza8s3zZ1lXL6zcyl++fZ3ALDPP0v7GBmoxUBRc4HKTt5yD6O2L2/HmWeiOjl6OR1h8ziC9qd/3qrlC7Bt4xrs2Hui9vm2+0rKmdIl4nYxDELP2/FeJ01YB8qgNuqfo+eF6MewYXUb9g4Mz1ifVtuJO/DaFRf9e57mWGbCxtL/sd5VgTluQaR9jAzUCAnAtdzDMDlqcX0ucYP9J0cCLYH9MkYb2aL/87zXe/qHceTMRRwduuScTr/oeNcAQGj75yBZYdxyQ5s26lbWe/qHazMm41eYj0RI0THlhWX1fbfJUYvStrSPkYEaIYQ4jheEmSQZuowRsJMt+j/P+z02MYmjQ5fiaTQJhXcNokj5gmSFccsNbdqoW1kD00Y1azuWtLRvQkjOcCkHUQKW4/zchGCgRkgAruUeJlUU27XjJDNpJDH0yyJtZIv+z/NeHxocqUlFNvW0x9R6YkMrcrAgWWEWckNdKqr3p0fW3YqrU5XEffYtQtwnimmGnmuYdbkMv/Tx+cOvYfUtC/HYXavq6q1t3z0wo71xH6OpyDtunG11w8VAjZSCKHlYrudq+dsXNdcsynHq+5p9S/d7Qn8AiQ1dxgiEl80B0/9Ezl283FQSknYtONIcXW6oX/ug9Umi39zpcseocqNM3TAJKTF6aRLb0htZ1xzTMUkU7+xcVlfGwysxAyR7jP5yIgBw403LVtgcBwM1UgqKkIfV7BjSPEb98+WGGzmOZIhJFhnVVdJGnhbVgZIkhy439PcD0/ok0eWWdf0pokQoCzdMQogPl2SMUbBpfxrHGGEfVjdYInIawDiAtwFMKaV6RWQ5gG8B6AJwGsBjSqmL4ZtACCEkKn4ZYyuukpt7VjSVp0V1oCTJEeRMmZZjpY4ut9TljrrTWhjpYxZumISQmdLBF/sv5Eq2bNP+tI5xU087jp0fw+aeFbVx8MyVn1nFTGGehP+qUuoN7fWTAF5SSn1ORJ6svv50iM8jJDWKkIfV7BjSPEZ9X+r621NNNi8NeZUF6m5+Jsc+XYffyIHSJSiZywZdbvli/wXjsr9vmZY9l8pW3DAJIXYYc6iAuu9v3r5//rHI1H7TNjv2nsC2jWsajlHjVyaxaN5s67w20//YG+YuWGxzHK1Ilh4G8N7q8tcA/AAM1IijuJ5vZkOzY0jzGPV9yec/+KPUduw4eZUF6m5++j9pU1HPRg6ULkHJXDbocsugIAww9y192XOpbMUNkxBih2m83HJ/F7bc3wUge2OQKPjHIlP7/dvs2HuiFqyZxiV92f/axmlZ/x/72S9ceiPwDRq2gZoCsFdEFID/qJR6FsAKpdT56t+HAVglxRGSBiziXMF/HgDMOC88V/GRV1mgyc1Pd+zTZWd5KXJNyVw2NEq0D+pb3o2Svmwqpp2XJ/mE5A19vNQDl7DFoF3CxvTDv822jWuwY++JwHHJW758bQoL5syy/v9iGsc++/aklZmAbaD2gFLqnIjcAmCfiPyz/kellKoGcTMQka0AtgJAZ2en5e4IaY0imIfEgU2RbJ6r+HBZFhhWlqk7Zj1/+DUcHbpUO648SDwpmYtO0rLRoEK0z///7Z15lBXVmcB/H20rO6IsioKNQdwjKKARjTJmyJEsZiaKOiYRzCQziWMco0mcMUbIyaZZTFwy0SSIGteoKOMWDYaoRAURFBBRhk2CiCKrCEL3N3/cW923i3qv6+3vdX+/c+r0fVV3r6rbde/9lpfWZLVSaeKshlEaHluwlheWzOw8UgAAGKZJREFUv8fmDzq2Y/qkfsjkLDv8XXGH16r6d/93nYhMA0YBb4vI/qr6lojsD6zLkPZm4GaAESNGJE7mDMMw2gPVLBZYiFjmkP7dOeXQfs3tqgURTxOZy5+yio0GVtCG9O/O8EG9gezO3U2c1TBKRK1bdywWmRxkx/unGhxei0g3oJOqbvHhscD3genA+cBP/N+HSllRw8iF9mA8pBikcZJtfVU8qlksMFexzNAa1vgRA1utFtaCiKeJzOVPqcVG45bWIito8ees3PUyjI7KBScNpmeX+lY6VKFj51KSj2PtUuWTZJ0xHKMix9n79+rMfr26lEV3L82OWn9gmohE8e9U1cdFZA5wr4h8GVgJjC9ZLQ0jRyqpZ1UMna9y6o2ZTlr7JC6emKtYZiie9tiCta3+6SXlVQvikIYj/KCJrCuGlFpsNC5Wm2QNMt96heKRdT36DChapQ2jHVPJha18HGuXKp+0oo/DB/VuzrfU/dXmRE1VlwHHJJxfD5xWikoZRi1TDJ2vYumNpdFRM9oncfHEfMUyI6fFbeVVC+KQhiN0Sh1ZVwwpp9hoKO4IFFyvUDyyrtve+5es4oZhFJ9iiRIWmk9a0ccyUIh5fsMwDKNKiYsn5iqWGYmAXDrWpXn1rc1Z86oFcUjDETqlTrpfpV5dzyRWG1qDzLdeoXjkAzu3by165Q3DKCrFcjpdjHzSiD6W2/G3TdQMo8gUQ+erWHpjaXTUjPZFJPq1+YOdBYmtheKNQKIj0FA8rZotXhqtxR2BrNYVS03oZDYsPzwfiUGGjmWTHLBnE4/sVN+5e7naZBhGfqRxTF2MfJKcesdFrJNEH0Px7CjfHz2ymOGDeqfWhUt0KF5Xn+ojzCZqhlFkiqHzVSy9MdM/63hEol8TT2zg4tMOyVtsLS7eGFnaCx2BhuJp1Wzx0mgt7hg6sq3E/QqdzIblh+ehRQwSMjudzSYeef+3Nr6FYRhVTRrH1MXIJ8lqbJKINZDRGm3oFHv+6k2pHV4nlV3XtVefNO2yiZph1ChpDI4UK062NPX9Bu+mw2pUjkj0K+6sNFfi4o1JjkBD8bRqtnhptBZ3LPTZKJS4k9mk85EYZOhYNnQ6C22LRzZueXdNqdpgGEZxSOOYuhj5JDn1jo8hSaKPcWu0l44dylUPLQKyj0Vtlb1qx7bNadrVKU0kwzCqjzQGR4oVJ1sa6VRX0gUfEZkiIutEZGHs/EUi8pqILBKRa4Lz/yUiS0VkiYh8spR1q0biIovZePO9bVz75Ou8+d62rOegRSQksnq1eO2WVGUY1UEkFlSs+xY9J7OXr098XpLiT56+iMnTF7UZF1osQ6587wMWr93S6rkLxZoqjY1PhlEbJP0PaxZDzDFOdO3u2at5Yfl7u8VJU3anzt17Z03ksR01wzCqnanADcBt0QkRGQOcARyjqjtEpJ8/fwRwDnAkMAD4s4gMVdXGste6QuQighiZMw5NGCedA9iy3a0Cvv/hLgAO36+HiTnWEJnEDfMlek6eWLQ2lRns0HQ20KbJ7Oh5iwifu0vHDq0mJ+ZTsfHJMIpGGl2ygkhjuTGtdUcfb96qDcxevr6VHm0YjoycROMYgHTqVJemCJuoGUaNksbgSLHiZEujTY272oheEKr6tIg0xE5/DfiJqu7wcdb582cAd/vzy0VkKTAKeK6UdawmSiWC2KNzPQDd9nL/Nk44eF/zl1ZDZBI3LJgSmauOnreI8LkbNXjfqjFYY+OTYRSXnHTJciCNU++0cSJCnbVQfzuuUxsZI4nGMQBtakq1QGMTNcMoMaVyXp0mn3zixOsLlM35dg4MBU4WkR8C24HLVHUOcADwfBBvtT9nJBD9Uwp3JsJzoRPrSHb/ix8bxNwVGwG38hlN1szhdcci/KBJs7sVftyE4UzOt6NV6J2NTQwbuDfjRw5s/nCavXx9RktrVeLw2sYnw8iTNLpk+ZDGvUfaOOGiV5L+dqZwaOp/1fatG9LU2yZqhlFiiuW8ulykcZJdTh21DOwB7AOcAIwE7hWRg3PJQES+CnwVYNCgQUWvYC2QtPsWnps8fVGzGCQ4E/1H7N8TaLHI15bIpNE+ydXXWqYdvUwikZEeCMDJh/Rl1GC3wn3L31bw/LL1GcUtQ5HMCjq8tvHJMPIkyUT+YwvWVs0uepyk+ka6aNnCaXXUzJiIYRi1yGrgAXXMBpqAPsDfgYFBvAP9ud1Q1ZtVdYSqjujbt2/JK2wYRpEopo5JabDxyTAKpbLvcO5IjuGU2I6aYZSYYjmvLhdpnGSXU0ctAw8CY4C/iMhQYE/gXWA6cKeI/AKnrH8IMLsC9asZsoksxkUj4+HRQ/Zl8nRnpvi4hr159a19OP3o/crbAKPihOKGaZy/hmQSiUwKJ+mPxBX4t2zfybADe/FPxw7gzzu3by24cflh41ONUHLDFUabhOLP0f+PiSc2tDLIEdcTqyaSTPqHIo7njDyQaS+tYUi/7q1EuCf9bMPbafK3iVo7o1T6UEZ6OsI9CNskV3/65VKWJSJ3AacCfURkNXAVMAWY4k1ifwicr6oKLBKRe4FXgV3AhWZRLTvZRBbjopHxcCQaCfD8sh7NpoyrVUTFKA25WoAMSeNTLTwXPYNxMUigVXjHriY61XfuXki70mDjU20Tf3YjTHy7fITiz9E7PPHEhqoyGpSNtKKPwwf1btWmSY07U+nC2EStnVFr+lDtkVq/B2l01MqJqp6b4dIXMsT/IfDD0tXISKTWxFSM4lOJZ6CIIkb5YONTO8HGr+qglu9DicYim6gZhmFUgGqwkhg5Hp54YkNeoiWhf5gLx3yEN95+v1p8W5WcTNYKOwqhyBi0iCqVy79ZJnGjUMToySs3vlXyihg1TdJzZOLb5SX8P3LOyAN5dMHbNXUP0og+FtImm6i1M2pNH6o9Uuv3II2OmlE41WAlMarDxBMb8ppohJb55q7Y2KHEhXJ14NzeiIuMlVtUKZu4USRi1Ljl3TVlqYxRs8xa2uKrK3qOZi1dXxMid+2FWUtbRJYfXfA2Lyx/r2z3IK4fl0bHNp7mR48sZv7qTa3GomkvrdmtTT9/4nUuHTu0WUeNuvpUH1Q2UWtntEd9qFqj1u9BrdffMIwyUmlRpQqIOxrth2j3NzQm0lGkAqqFM487sNkFTDl35cE5144W3EIn1dmca8fTzF+9yV0Ixp8h/bszfJCzvp/JKXZd11590tTRJmqG0U7I14hJPg6uwzT1/QYfU3DlOyBJzqartQ5xMc3o93ENezNvVS+G9OvO6Ufvx7VPvt5hHF5nslDYUcjV4XUa2hInDcUtQ+uOSRbVZi9fX0mH10aNkKs/QKP4xA0HlfM+jB6yL08scmPGuKP7N5/LJU2SM+vxIwa2akeSU+xVO7ZtTlNHm6gZRjshXyMm+RgPqQKH1zVPkrPpaq1DXEwz+n34fj2aRc0iK3wdxeF1JmuFHYVSfOC2JU6aZKFvx66mVuKOoTXICjq8NgyjBgjF9yPR6bYsFyeliYfjeSSJapvDa8MwDKO0mKiZUQnSiDvas2kYRi7kM2bkKnqdRxniXHuUBxF5B1gJ9ME5f6wVrL6lxepbBOr7DT4m2t3SpsZdO9ctj/ybZa1vPB207JLF8klMs2vjWm38YEtNL/r4sel9yn9fK/Us5VZuXf2edV179WnctuldGnd+GP1u2rFtc6e9uvZs3LbpXXAy981xCi2zONRG/1qZjrr6Peu69e4P0Pj+hrd3e47q6ves69JzUNOHH6yNVqObtm/d0PwMxp5NmhoHNu3c3q1o9asQIrIFWFKBomv/mareMitVbkcpM125wZiz21iSNk04FnXuNgDkg93Gr/B/po/fuOWdem1q2qetRpR1otZcqMiLqjqi7AXnidW3tFh9S0ut1bdSVKKfKnVvOkpbrX+tzPaAPcftr8xKldtRyqxUuaUos6ZXwQ3DMAzDMAzDMNojNlEzDMMwDMMwDMOoMio1Ubu5QuXmi9W3tNwMII6nRKRnIZn5fK4TkaUi8oqIHOvP9xWRx4tV3xqi1upbKSrRT5W6NwWVm+e7uluZInKYiDwnIjtE5LKUZc8UkbSiJSXrXxHZWq5yRWSCiNyQ4dqjIrJ3tjJz7LNcuFlEJqW9d0F9Pi0i38+3zDzT1TI1OU5YmVVZbkcps1LlFr3MiuioGe0TEdlDVXcVkP5TwCdU9ZIC6zEOuAgYBxwP/EpVj/fXbgF+p6qzCinDMGqZKnpX+wEHAZ8DNqjqz1KkmQlcpqovFlJ2LiT1l4hsVdXuKdML7v9tU4brK1S1IUv6CcAIVf2P9LVulX4mReizDP0wCdia5t5FeQCNwEvAaFXdVkidDKMY+Hd0BvA5VU3l3ypDPocBtwDHAlekfS8y5DUVeFhV78vlHU4zXohIX+BhYE/gG6r6TMo6TQCeUNU1aeKnRUSGAQNU9VH/exI5jCsZ8rwbuFJV3yhOLSuDiT7WKCLSICKLReS3IrJIRJ4QkS7+2jARed7vJk0Tkd7+/EwRuVpEZovI6yJyckK+A0RkfnA0ishBfjfqfhGZ44/RPv4kEbldRGYBt4tIZxG5RUQWiMg8ERnj4x3py53v63VIQrPOAx4K2veaiNzh23mfiKT1onsGcJs6ngf2FpHIn86DvhzDKAv2rmZGVdep6hxgZ559e66v/0IRudqfO0tEfuHDF4vIMh8+2Lc9nke2e/BLEXkRuFhEBovb/VsgIj+I5fEt39eviMjkoF+WiMhtwEJgYD5tDBggIo+LyBsick1Q9goR6ePDV/oynxWRu6T1TtdZbTxPIiI/9X25QETO9udPFZFnRGQ68Ko/d4XP51ng0CCPj/g6zvVpDvPnp4rIb0TkBeAadSvEM4FPF9gnhgE0LwAUwjjg5UImaZ73gG8AeU8wysRpwAJVHZ52kuaZAJTCkfww3D0oCiJSB/wP8O1i5VkxVNWOGjyABmAXMMz/vhf4gg+/Apziw98HfunDM4Gf+/A44M9tlHEhcK8P3wmc5MODgMU+PAmYC3Txvy8FpvjwYcAqoDNwPXCeP79nFD9W3kqgR9A+xa24AkzBrSYBXAvMTzgu99cfjurqf8/ArS4BHIAbnCp+D+3oGIe9q5nf1SC/SVGaFP05ExiB+1hYBfQF9gCewu3M7QfM8XHvA+b49/584McJ+WW7B78O4k0HvhT091YfHosTdxHc4ufDwMd9vzQBJ6Ro04o2rk8AlgG9/D1aCQyM0uLMUI/0fdsZ6AG8EdyHNp8n4PPAk0Ad0N/37f7AqTjXFYN9vOOABUBXoCewNChnBnCIDx8PPOXDU32/1AXlnQdcX+n3047cD/9sLwZ+CywCnqBlXBkGPO/fq2lA7+AZvBqYDbwOnJyQ7wBajxONuB33vsD9/l2eQ8tYMwm4HZgF3OWf/Vv88zkPGOPjHenLne/rdUhC2XcCpwbtew24w7fzPqBrjn00ifRj2vd8uxbixxJ/fipwZtB/I7LkMdH362x/X24I2vKUb/cM3P+EYf79fsf3STdf1kLfd5dkKONMYCvOVcR8oAtuwjfPp5sC7OXjjvN9OBe4DrcziC9riq/nPNzC+p6x+pzt+2+Kb/cy3K5fVI8vBPfzJvy44uv2c+Bl4CTceLwc2KPS70whh+2o1TbLVXW+D88FGkSkF7C3qv7Vn78V99EQ8UAYP1PGfhX+K8AF/tQngBtEZD7ug6WniERiP9NV9QMfPgn4A4Cqvob7oBgKPAf8t4h8BzgoiB+yj6puCX6/qS0iin/weaOql6jqsITjJ5naE7CO0qwGGUY27F3N/V1ti5HATFV9R5043h3Ax1V1LdBdRHrgdrHuxPXryUCrleMU9+CeIDwa9zEI7uMwYqw/5uHE+Q4Dol3Ilep29XdDRG6MdkNxu2XRzugVGdo7Q1U3qep23M7WQbHro4GHVHW7vzf/G7ve1vN0EnCXqjaq6tvAX3F9DDBbVZf78MnANFXdpm73YbpvT3fgROCPvk034SZ6EX9U1cbgt43Ftc0hwI2qeiSwETfRB7gN+I6qfhT38X5VkGYPVR0F/GfsPACquiYaI3CTjftVdSXwK+BaVR3py/ldkOwInBj2ubgFFFXVo4FzgVtFpDPw7zgViGG4RZ7VCe0ZjXs3Ig7FLdQcDmwGvg4gItfGJBmi4/J03ZbIDao6UlWPwk1+ctppFicxNNm34SRcn0RcD9zq78cdwHX+f9H3gHt8nxwGHKCqR/m+uyWpHFW9D3gRt5A3DLdANxU426fbA/ia7/ObgNNV9TjcRDviCtwCzihgDPBToD6sj6pG4+5hwCeBUcBVIlIvIofjJnKjfR0aaZGS6ga8oKrHqOqz6kTNlwLH5NKf1UahW8VGZdkRhBtxL3jaNI34+y9Ob2s4sEZVx/mX/vfAZ1U1UpjvhFsZ3h5mJiLgVluzoqp3erGXTwGPisi/qepTsWi7RKSTtuhxxBUo1Zd5Le4Fj3O3/wD8O63FjA7058CtuCV9eBpGKbF3tTXRu1oq/oZbYV6Cm5xdAHwMt4uYC/H+SlLqFtxO3U2tToo0JKRvyUj1wiDuCv/RkY34M5Tr/+/dnqccaPO5wT13G7O0I56HjcW1TdrFpz8GaXJdfDrJn/oEcIQfwyD74tP14BafRCRcfLpCRA4EHtBknaW2Fp++AfxMC9TLzcAYEfk2bpd6H9wuZXyhJRvH4xetAETkHly7wY17/+zDtwPX7J6cZcDBInI98AhuhzQNh+Keg9f971txk+WZwLJgcecu4Ks+PBb4bCCW3Rm3y5fEI6q6A9ghIutwO/2n4Xb15/jnoQtu0Qfc2HZ/LI9oQWguNYrtqLUzVHUTsCHQQfgibmU0W5qJfhVjnIjU4wbW7wQvH7gX96LohzjFzySewa9uiMhQ3Au4REQOxr241+F0Wz6akHYJcHDwe5CIfMyH/wV41te3rVX66cCXvM7FCcAmVX3LXxuK2943jIpi72p2RGSGiByQJcps4BQR6eP1Ec6lpf+eAS4DnsaLQAE7fJ83k+M9mAWc48OhnuufgAuiD0cROUCckZRyMwv4jDjdw+7krv/1DHC2iNSJMzTwcVwfx3ka+JyIdPG7lp8B8Ltry0XkLGjWecu2km1jcW2Tz8JB4uKT35GKjEhEi0/jExafovHjgOBaqsUn4LO4hYFHReQfEqLtEpHwmzjj4lMxd9T87tOvcSKOR+N2Ejvnk1e+qOoG3K7TTNzu4++yJigMAT4f3MtBqro4Q9ykZ0xwO4RR+kNVdZKPsz22aw/tYEHIJmrtk/OBn4rIKzhZ5FzMIJ+IEw2YHAxAA3CrSSPEKcu/inuZk/g10ElEFuDEhib4FZHxwEIvEnMUTjwiziM4fYiIJcCFIrIY6I1TDE3Do7gVoqW4Qe/rwbUxvhzDqAY69LsqIvuJyGrgm8B3RWS1iPT0H0xDcIr5ifjFl8uBv+B0Euaq6kP+8jO4XfWn/T/uN/GTxwTS3oOLfRsX4HTeono8gROvfM5fuw+nI1ZW1BllmY7TRXkMJ3a2KWui1kzzaV/G6bR8W50Yabycl3DPy8u+nDnB5fOAL4vIy7hdgTOylGdjcTvDFp+yk2HxKZqUvesXWM5sK58EXsAtWu3r+/Cs4NrfaL3AtJvhEHHGiDqp6v3Ad3EWKzOxhZbxbQluJ3WI/x3d7yW4HboGf/7sIP2fgIvEb4eJyPCEfLMxAzgzWgwTkX1EJC4GHlL7C0JaBYpydtihquD0GZ704QZgYQnKeBqv3GyHHXbkd5T6XcVNEH9R6XbW2gF093+74nRJjq10nTLUsz9O567idbEjr/vX6p3H7V5P8uHQmMiDtDYmEhn16kOCAR3gFGA7rQ2KDPDx7/F5vgr8xsefRGCwg8zGRC7HLRzMBx7HiTnGy74S+Negfa/hRB4X48TpUhkTwRkzWo3Ta9vowz1xGyMrSTbO9APg/3C74rcEfTmV/IyJ3EyLMZGDiBkT8ecnBHGOwenXRn1+epZyPk86YyKfocWYyG+AO/z5Ljj9tQX+nkRGRvbBLfqExkTCe7sQaPDhs2kxDDMXb7AJb+ApSNMfp19b8XemkMP8qBlVhYiMxw+kuBf4qCLm3RengPpgsfI0jI5KKd9VIz9E5E6cIYHOOPGgH1e4SomIyEhgp7boOBlGRfEil7ep6j/6naBif38cBVygqt8sVp7VjIh0V9WtfufsRuANVb22zHW4BNisqr8vZ7nFxiZqhmEYhmEYRofGFp+Kh58knY8zvT8P+IqW2bm9iEwEbldnFbhmsYmaYRiGYRiGYVQ54izy7hU7/UVVXVDkcm7EmfsP+ZWqJpruN0qHTdQMwzAMwzAMwzCqDLP6aBiGYRiGYRiGUWXYRM0wDMMwDMMwDKPKsImaYRiGYRiGYRhGlWETNcMwDMMwDMMwjCrDJmqGYRiGYRiGYRhVxv8DBYW6TMauFV8AAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 1080x1080 with 3 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "try:\n",
    "    V0 = L2(mesh,order=0,dgjumps=True)\n",
    "    a0 = BilinearForm(V0)\n",
    "    a0.Assemble()\n",
    "\n",
    "    V1 = L2(mesh,order=1,dgjumps=True)\n",
    "    a1 = BilinearForm(V1)\n",
    "    a1.Assemble()\n",
    "\n",
    "    V2 = L2(mesh,order=1,dgjumps=True, all_dofs_together=True)\n",
    "    a2 = BilinearForm(V2)\n",
    "    a2.Assemble()\n",
    "\n",
    "    import scipy.sparse as sp\n",
    "    import matplotlib.pylab as plt\n",
    "    plt.rcParams['figure.figsize'] = (15, 15)\n",
    "\n",
    "    a0.mat.AsVector()[:] = 1 # set every entry to 1\n",
    "    rows,cols,vals = a0.mat.COO()\n",
    "    A0 = sp.csr_matrix((vals,(rows,cols)))\n",
    "    a1.mat.AsVector()[:] = 1 # set every entry to 1\n",
    "    rows,cols,vals = a1.mat.COO()\n",
    "    A1 = sp.csr_matrix((vals,(rows,cols)))\n",
    "    a2.mat.AsVector()[:] = 1 # set every entry to 1\n",
    "    rows,cols,vals = a2.mat.COO()\n",
    "    A2 = sp.csr_matrix((vals,(rows,cols)))\n",
    "    fig = plt.figure()\n",
    "    ax1 = fig.add_subplot(131)\n",
    "    ax2 = fig.add_subplot(132)\n",
    "    ax3 = fig.add_subplot(133)\n",
    "    ax1.set_xlabel(\"non-zeros (p=0)\")\n",
    "    ax1.spy(A0,markersize=3)\n",
    "    ax2.set_xlabel(\"non-zeros (p=1, low order + high order)\")\n",
    "    ax2.spy(A1,markersize=1)\n",
    "    ax3.set_xlabel(\"non-zeros (p=1, all_dofs_together)\")\n",
    "    ax3.spy(A2,markersize=1)\n",
    "except ImportError:\n",
    "    pass"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "outputs": [],
   "source": [
    "try:\n",
    "    plt.show()\n",
    "except ImportError:\n",
    "    pass"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## Supplementary 1: Skeleton formulation\n",
    "\n",
    "So far we considered the DG formulation with integrals on the boundary of each element. Instead one could formulate the problem in terms of facet integrals where every facet appears only once. \n",
    "\n",
    "Then, the corresponding formulation is:\n",
    "\n",
    "Find $u: [0,T] \\to V_h := \\bigoplus_{T\\in\\mathcal{T}_h} \\mathcal{P}^k(T)$ so that\n",
    "$$\n",
    "  \\sum_{T} \\int_T \\partial_t u v + b \\cdot \\nabla u v + \\sum_{F\\in\\mathcal{F}^{int}} \\int_{F} b_n (u^{neighbor} - u) v^{downwind} + \\\\\n",
    "  \\sum_{F\\in\\mathcal{F}^{inflow}} \\int_{F} b_n (u^{inflow}-u) v = \\int_{\\Omega} f v , \\quad \\forall v \\in V_h.\n",
    "$$\n",
    "\n",
    "Here $\\mathcal{F}^{int}$ is the set of interior facets and $\\mathcal{F}^{inflow}$ is the set of boundary facets where $b \\cdot \\mathbf{n} < 0$. $v^{downwind}$ is the function on the downwind side of the facet and $u^{inflow}$ is the inflow boundary condition. \n",
    "\n",
    "In NGSolve these integrals can be assembled with `skeleton=True`. The facet integrals are divided into interior facets and exterior (boundary) facet. Hence, we add `BND`/`VOL`:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "outputs": [],
   "source": [
    "b = CoefficientFunction((1,2))\n",
    "n = specialcf.normal(mesh.dim)\n",
    "ubnd = IfPos(x,1,0)\n",
    "\n",
    "V = L2(mesh,order=2)\n",
    "u,v = V.TrialFunction(), V.TestFunction()\n",
    "\n",
    "c = BilinearForm(V)\n",
    "c += SymbolicBFI( b * grad(u) * v)\n",
    "\n",
    "bn = b*n\n",
    "vin = IfPos(bn,v.Other(),v)\n",
    "c += SymbolicBFI( bn*(u.Other() - u) * vin, VOL, skeleton=True)\n",
    "c += SymbolicBFI( IfPos(bn, 0, bn) * (u.Other(ubnd) - u) * v, BND, skeleton=True)\n",
    "\n",
    "gfu_expl = GridFunction(V)\n",
    "Draw(gfu_expl,mesh,\"u_explicit\")\n",
    "\n",
    "res = gfu_expl.vec.CreateVector()\n",
    "c.Apply(gfu_expl.vec,res)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "outputs": [],
   "source": [
    "t = 0\n",
    "dt = 0.001\n",
    "tend = 1\n",
    "\n",
    "while t < tend-0.5*dt:\n",
    "    c.Apply(gfu_expl.vec,res)\n",
    "    V.SolveM(res,rho=CoefficientFunction(1.0))\n",
    "    gfu_expl.vec.data -= dt * res\n",
    "    t += dt\n",
    "    Redraw(blocking=True)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## Supplementary 2: DG for Diffusion  \n",
    "\n",
    "DG formulation for the Poisson problem\n",
    "$$\\text{find: } u \\in H_{0,D}^1 \\quad \\int_\\Omega \\nabla u \\nabla v = \\int_\\Omega f v \\quad \\forall v \\in H_{0,D}^1$$\n",
    "\n",
    "For the non-conforming discretization we use the $L^2$ space and the (symmetric) interior penalty method:\n",
    "Find $u \\in V_h := \\bigoplus_{T\\in\\mathcal{T}_h} \\mathcal{P}^k(T)$ so that\n",
    "\\begin{align}\n",
    "  & \\sum_{T} \\int_T \\nabla u \\nabla v + \\sum_{F\\in\\mathcal{F}^{int}} \\int_{F} \\{\\!\\!\\{ - \\nabla u \\cdot n \\}\\!\\!\\} [\\![v]\\!] +  \\{\\!\\!\\{ - \\nabla v \\cdot n \\}\\!\\!\\} [\\![u]\\!] +  \\frac{\\alpha}{h} [\\![u]\\!][\\![v]\\!] \\\\\n",
    "  & + \\sum_{F\\in\\mathcal{F}^{ext}} \\int_{F} - \\nabla u \\cdot n v - \\nabla v \\cdot n u + \\frac{\\alpha}{h}~ u v = \\sum_{F\\in\\mathcal{F}^{ext}} \\int_{F} - \\nabla v \\cdot n u_{D} + \\frac{\\alpha}{h} u_{D} v  + \\int_{\\Omega} f v , \\quad \\forall v \\in V_h.\n",
    "\\end{align}\n",
    "where $\\{\\!\\!\\{\\cdot\\}\\!\\!\\}$ and $[\\![ \\cdot ]\\!]$ are the usual average and jump operators across facets."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "outputs": [],
   "source": [
    "n = specialcf.normal(mesh.dim)\n",
    "h = specialcf.mesh_size\n",
    "IfPos(x-y,0,1)\n",
    "source = IfPos(x-y,5,-5)\n",
    "\n",
    "order=2\n",
    "V = L2(mesh,order=order,dgjumps=True)\n",
    "a = BilinearForm(V)\n",
    "u,v = V.TrialFunction(), V.TestFunction()\n",
    "a += SymbolicBFI( grad(u) * grad(v))\n",
    "def avg_flux(u):\n",
    "    return 0.5*(-grad(u)*n-grad(u.Other())*n)\n",
    "def jump(u):\n",
    "    return u-u.Other()\n",
    "alpha = 5 * order * (order+1)\n",
    "a += SymbolicBFI( avg_flux(u) * jump(v) + avg_flux(v) * jump(u) + alpha/h * jump(u) * jump(v), VOL, skeleton=True)\n",
    "a += SymbolicBFI( - grad(u)*n*v - grad(v)*n*u + alpha/h * u * v, BND, skeleton=True)\n",
    "a.Assemble()\n",
    "\n",
    "f = LinearForm(V)\n",
    "f += SymbolicLFI( (- grad(v)*n + alpha/h * v) * ubnd, BND, skeleton=True)\n",
    "f += SymbolicLFI( source * v)\n",
    "f.Assemble()\n",
    "\n",
    "gfu_impl = GridFunction(V)\n",
    "gfu_impl.vec.data = a.mat.Inverse() * f.vec\n",
    "Draw(gfu_impl,mesh,\"u_implicit\")"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.6.5"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}
